{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pima Indians Diabetes Data Set - Logistic回归分析\n",
    "\n",
    "数据说明：\n",
    "Pima Indians Diabetes Data Set（皮马印第安人糖尿病数据集） 根据现有的医疗信息预测5年内皮马印第安人糖尿病发作的概率。   \n",
    "\n",
    "数据集共9个字段: \n",
    "0列为怀孕次数；\n",
    "1列为口服葡萄糖耐量试验中2小时后的血浆葡萄糖浓度；\n",
    "2列为舒张压（单位:mm Hg）\n",
    "3列为三头肌皮褶厚度（单位：mm）\n",
    "4列为餐后血清胰岛素（单位:mm）\n",
    "5列为体重指数（体重（公斤）/ 身高（米）^2）\n",
    "6列为糖尿病家系作用\n",
    "7列为年龄\n",
    "8列为分类变量（0或1）\n",
    "\n",
    "数据链接：https://archive.ics.uci.edu/ml/datasets/Pima+Indians+Diabetes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import必要的工具包，用于文件读取／特征编码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  # 矩阵操作\n",
    "import pandas as pd # SQL数据处理\n",
    "\n",
    "from sklearn.metrics import r2_score  #评价回归预测模型的性能\n",
    "\n",
    "import matplotlib.pyplot as plt   #画图\n",
    "import seaborn as sns\n",
    "\n",
    "# 图形出现在Notebook里而不是新窗口\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "数据文件路径和文件名"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pregnants</th>\n",
       "      <th>Plasma_glucose_concentration</th>\n",
       "      <th>blood_pressure</th>\n",
       "      <th>Triceps_skin_fold_thickness</th>\n",
       "      <th>serum_insulin</th>\n",
       "      <th>BMI</th>\n",
       "      <th>Diabetes_pedigree_function</th>\n",
       "      <th>Age</th>\n",
       "      <th>Target</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>0.639947</td>\n",
       "      <td>0.866045</td>\n",
       "      <td>-0.031990</td>\n",
       "      <td>0.670643</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>0.166619</td>\n",
       "      <td>0.468492</td>\n",
       "      <td>1.425995</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>-0.844885</td>\n",
       "      <td>-1.205066</td>\n",
       "      <td>-0.528319</td>\n",
       "      <td>-0.012301</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>-0.852200</td>\n",
       "      <td>-0.365061</td>\n",
       "      <td>-0.190672</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>1.233880</td>\n",
       "      <td>2.016662</td>\n",
       "      <td>-0.693761</td>\n",
       "      <td>-0.012301</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>-1.332500</td>\n",
       "      <td>0.604397</td>\n",
       "      <td>-0.105584</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>-0.844885</td>\n",
       "      <td>-1.073567</td>\n",
       "      <td>-0.528319</td>\n",
       "      <td>-0.695245</td>\n",
       "      <td>-0.540642</td>\n",
       "      <td>-0.633881</td>\n",
       "      <td>-0.920763</td>\n",
       "      <td>-1.041549</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>-1.141852</td>\n",
       "      <td>0.504422</td>\n",
       "      <td>-2.679076</td>\n",
       "      <td>0.670643</td>\n",
       "      <td>0.316566</td>\n",
       "      <td>1.549303</td>\n",
       "      <td>5.484909</td>\n",
       "      <td>-0.020496</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pregnants  Plasma_glucose_concentration  blood_pressure  \\\n",
       "0   0.639947                      0.866045       -0.031990   \n",
       "1  -0.844885                     -1.205066       -0.528319   \n",
       "2   1.233880                      2.016662       -0.693761   \n",
       "3  -0.844885                     -1.073567       -0.528319   \n",
       "4  -1.141852                      0.504422       -2.679076   \n",
       "\n",
       "   Triceps_skin_fold_thickness  serum_insulin       BMI  \\\n",
       "0                     0.670643      -0.181541  0.166619   \n",
       "1                    -0.012301      -0.181541 -0.852200   \n",
       "2                    -0.012301      -0.181541 -1.332500   \n",
       "3                    -0.695245      -0.540642 -0.633881   \n",
       "4                     0.670643       0.316566  1.549303   \n",
       "\n",
       "   Diabetes_pedigree_function       Age  Target  \n",
       "0                    0.468492  1.425995       1  \n",
       "1                   -0.365061 -0.190672       0  \n",
       "2                    0.604397 -0.105584       1  \n",
       "3                   -0.920763 -1.041549       0  \n",
       "4                    5.484909 -0.020496       1  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#input data\n",
    "fileName = r\"D:\\AllFile\\AI\\jupyter_notebook\\Logistic_Regression\\FE_pima-indians-diabetes.csv\"\n",
    "train = pd.read_csv(fileName)\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  数据基本信息\n",
    "样本数目、特征维数\n",
    "每个特征的类型、空值样本的数目、数据类型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 768 entries, 0 to 767\n",
      "Data columns (total 9 columns):\n",
      "pregnants                       768 non-null float64\n",
      "Plasma_glucose_concentration    768 non-null float64\n",
      "blood_pressure                  768 non-null float64\n",
      "Triceps_skin_fold_thickness     768 non-null float64\n",
      "serum_insulin                   768 non-null float64\n",
      "BMI                             768 non-null float64\n",
      "Diabetes_pedigree_function      768 non-null float64\n",
      "Age                             768 non-null float64\n",
      "Target                          768 non-null int64\n",
      "dtypes: float64(8), int64(1)\n",
      "memory usage: 54.1 KB\n"
     ]
    }
   ],
   "source": [
    "train.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 从原始数据中分离输入特征x和输出y\n",
    "y_train = train['Target']   \n",
    "X_train = train.drop([\"Target\"], axis=1)\n",
    "\n",
    "#保存特征名字以备后用（可视化）\n",
    "feat_names = X_train.columns \n",
    "\n",
    "# 对稀疏矩阵进行压缩处理：把稀疏矩阵按照某种格式处理，可以压缩数据量，比如只记住非零值的数值和位置，没有保存的就都是0。\n",
    "from scipy.sparse import csr_matrix\n",
    "X_train = csr_matrix(X_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 交叉验证用于评估模型性能和进行参数调优（模型选择）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 默认参数的Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logloss of each fold is:  [0.48797856 0.53011593 0.4562292  0.422546   0.48392885]\n",
      "cv logloss is: 0.47615970944434044\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "lr = LogisticRegression()\n",
    "# 忽略警告信息\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "#分类任务中交叉验证缺省是采用StratifiedKFold\n",
    "#数据集比较大，采用3折交叉验证\n",
    "from sklearn.model_selection import cross_val_score\n",
    "loss = cross_val_score(lr, X_train, y_train, cv=5, scoring='neg_log_loss')\n",
    "print ('logloss of each fold is: ',-loss)\n",
    "print ('cv logloss is:', -loss.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression + GridSearchCV\n",
    "\n",
    "logistic回归的需要调整超参数有：C（正则系数，一般在log域（取log后的值）均匀设置候选参数）和正则函数penalty（L2/L1） \n",
    "目标函数为：J =  C* sum(logloss(f(xi), yi)) + penalty \n",
    "\n",
    "在sklearn框架下，不同学习器的参数调整步骤相同：\n",
    "1. 设置参数搜索范围\n",
    "2. 生成学习器实例（参数设置）\n",
    "3. 生成GridSearchCV的实例（参数设置）\n",
    "4. 调用GridSearchCV的fit方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 采用5折交叉验证，用log似然损失，对Logistic回归模型的正则超参数调优"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "#需要调优的参数\n",
    "# 请尝试将L1正则和L2正则分开，并配合合适的优化求解算法（slover）\n",
    "# tuned_parameters = {'penalty':['l1','l2'],\n",
    "#                   'C': [0.001, 0.01, 0.1, 1, 10, 100, 1000]\n",
    "#                   }\n",
    "penaltys = ['l1','l2']\n",
    "Cs = [0.001, 0.01, 0.1, 1, 10, 100, 1000]\n",
    "tuned_parameters = dict(penalty = penaltys, C = Cs)\n",
    "\n",
    "lr_penalty= LogisticRegression(solver='liblinear')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, error_score='raise-deprecating',\n",
       "             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,\n",
       "                                          fit_intercept=True,\n",
       "                                          intercept_scaling=1, l1_ratio=None,\n",
       "                                          max_iter=100, multi_class='warn',\n",
       "                                          n_jobs=None, penalty='l2',\n",
       "                                          random_state=None, solver='liblinear',\n",
       "                                          tol=0.0001, verbose=0,\n",
       "                                          warm_start=False),\n",
       "             iid='warn', n_jobs=4,\n",
       "             param_grid={'C': [0.001, 0.01, 0.1, 1, 10, 100, 1000],\n",
       "                         'penalty': ['l1', 'l2']},\n",
       "             pre_dispatch='2*n_jobs', refit=True, return_train_score='true',\n",
       "             scoring='neg_log_loss', verbose=0)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 如果 return_train_score=\"false\"，cv_results_属性将不包括训练分数。\n",
    "grid= GridSearchCV(lr_penalty, tuned_parameters,cv=5, scoring='neg_log_loss',n_jobs = 4, return_train_score=\"true\")\n",
    "grid.fit(X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4760264922287083\n",
      "{'C': 1, 'penalty': 'l1'}\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print(-grid.best_score_)\n",
    "print(grid.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUVf7/8ddnJr1QEzD0rtRQsoKg2BARIXQERYoF1wboLtb9qaCrsCrifsUVQRFRKRaKiDQRUKQFDL0jJbTQQklP5vz+mAmEMEDacDPJ5+ljHpm5c+6d9wWcT24554gxBqWUUionm9UBlFJKFU1aIJRSSrmlBUIppZRbWiCUUkq5pQVCKaWUWz5WBygsYWFhpkaNGlbHUEopr7Ju3boTxphwd+8VmwJRo0YNYmJirI6hlFJeRUT2X+k9j55iEpEOIrJDRHaLyEtu3v9ARGJdj50ikpDtvQEissv1GODJnEoppS7nsSMIEbED44B7gDhgrYjMMcZszWpjjHkuW/tngWau5+WA14EowADrXOue9lRepZRSl/LkEcTNwG5jzF5jTBowDehylfZ9gamu5/cCi4wxp1xFYRHQwYNZlVJK5eDJaxCVgYPZXscBLd01FJHqQE1gyVXWrexmvcHAYIBq1aoVPLFSyiulp6cTFxdHSkqK1VGKrICAAKpUqYKvr2+u1/FkgRA3y6408FMf4DtjTGZe1jXGfAp8ChAVFaWDSilVQsXFxREaGkqNGjUQcff1UbIZYzh58iRxcXHUrFkz1+t58hRTHFA12+sqwOErtO3DxdNLeV1XKVXCpaSkUL58eS0OVyAilC9fPs9HWJ4sEGuBuiJSU0T8cBaBOTkbiciNQFlgZbbFC4D2IlJWRMoC7V3LlFLKLS0OV5efPx+PFQhjTAbwDM4v9m3ADGPMFhEZKSLR2Zr2BaaZbOOOG2NOAW/iLDJrgZGuZYUuMyODVZ88xeF9OzyxeaVUEfXA+JU8MH7ltRuWYB7tB2GMmWeMqWeMqW2M+bdr2WvGmDnZ2rxhjLmsj4Qx5nNjTB3XY5KnMh7+awsNjs7C94t72bVxlac+RilVzIWEhFx43qFDB8qUKUOnTp3ctn366adp2rQpDRo0IDAwkKZNm9K0aVO+++67PH3m+vXrmT9/foFyX02JH4upat1ITvWejUGo+H03/vztR6sjKaW83PDhw5kyZcoV3x83bhyxsbHMmzeP2rVrExsbS2xsLD179szT52iBuA5qNPgb8tgiEuzlabB4ECvmfGZ1JKWUF7v77rsJDQ3N17q7du3i3nvvpUWLFrRt25adO3cCMG3aNBo1akRkZCR33nknycnJjBw5kq+//jpfRx+5UWzGYiqo8Cp1OD/kVw58HM0t6/7B4oRj3P3wy3rhSykvM+LHLWw9fPaa7bYecbbJzXWIBpVK8XrnhgXOlhuDBw9m4sSJ1K5dmxUrVvDMM8+wcOFCRowYwdKlS6lYsSIJCQkEBgby2muvsXnzZsaOHeuRLFogsgkpE06N5xexbVwv2u0dzaJxR7jj72Px9bFbHU0pVQIkJCSwatUqevTocWFZRkYGAG3atKF///706tWL7t27X5c8WiBy8A0IocGwOWwa/yj3xH/Jrx8cJ+rpSYQGBVodTSmVC7n9TT/ryGH6E7d4Mk6eGGMICwsjNjb2svcmTJjA6tWrmTt3LpGRkWzcuNHjefQahBti96Xxk5PZUmcwdyb+zOYPunDspEfuslVKqQvKli1LREQEM2fOBMDhcLBhwwYA9u7dS6tWrXjzzTcpW7Yshw4dIjQ0lHPnznksjxaIKxGhYb932Rn1Bi3T1nDsow7s3n/A6lRKKS9w22230atXL3755ReqVKnCggW57+c7bdo0PvnkEyIjI2nYsCFz584F4LnnnqNx48Y0btyYdu3a0ahRI+666y42bNhAs2bNPHKRWrL1T/NqUVFRxlMTBh34fSo3LH6Gg1Qkoft0WjRp7JHPUUrlz7Zt26hfv36e1imKp5g8zd2fk4isM8ZEuWuv1yByodqtfYkvFc4NPzxM8PfR/HpmEnfedofVsZRSBVCSCkN+6SmmXKrQpB2OgT/jZ4fmi/sye/a3FJejL6WUckcLRB6E1mhK0FO/kuxXng7rn2Tql/8j06FFQilVPGmByKOAsBpUGLqUEyH1eGDvK3zz8Rskp2Vee0WllPIyWiDywRYSRuWhizgc3oaHT4xl9thnOHlOZ7JSShUvWiDyyy+Yqk/OIq56d/okfcOKDx9mX/y1u/crpYqISfc7H+qKtEAUhN2XKgM/50iTp4jOWMjej3sQu/eI1amUUha43sN9z5w5k3fffbfAua9Gb3MtKBEiur/DidIR3PHba6yf3IUlXadwV7MbrU6mlLLI8OHDSUpKYvz48W7fHzduHAD79u2jU6dObofWAOc4TD4+7r+mu3XrVjhhr0KPIApJ2N1DON/pUyJlD5Vnduf7X1dbHUkpZZGCDPd966238uqrr9K2bVs++ugjZs+eTcuWLWnWrBnt27cnPj4egIkTJzJs2DAA+vXrx9ChQ2ndujW1atW6MFRHQekRRCEqFdWblFJhVJv6ECFL+zLh1Ec82u0+bDYdMlyp6+bnl+Dopmu3O+oa7C431yFuaAz3jSpYrjw4e/Ysy5cvB+D06dNER0cjInzyySe8//77jB49+rJ14uPjWbFiBZs2baJ3796FcoShRxCFLKDeXfg+Np9SvtBr42N8+MVXpGbobbBKqdzr06fPhecHDhygffv2NG7cmDFjxrBlyxa363Tt2hURoUmTJhw6dKhQcugRhAf4VI4k5KlfODMhmif3P8d/Pz7O4MefpXSgr9XRlCr+cvubftaRw6CfPJcln4KDgy88f/rpp3nllVfo2LEjixcvZtQo9/vn7+9/4XlhjfKgRxAeIuVqUuaZX0kqeyPPnxzJZ/99g0MJyVbHUkp5mTNnzlC5cmWMMUyePPm6frYWCE8KDqPckws4W+lWnk/+iB//7zm2HEqwOpVSysMKMtx3Tm+88QbdunXj9ttvp2LFioWY8tp0uO/rITOdM9OfoPTO75lq2lPlwf/jthtvsDqVUsVGfob7LsqnmDwlr8N96xHE9WD3pXSfiZxv8RR9ZSGJXz/MD2v2WJ1KqZJt0E8lqjjkhxaI68VmI6TzO6Tc9SYdbGuoNPchPlmwXocMV0oVWVogrrOAtkPI6DqBKNsubl/Rn7en/0p6psPqWEopdRktEBbwadobe79vqe1zkgHbHue1z34gMTXD6lhKKXUJLRAWkTp34ffYPMr7Oxh+aCivjZtEvA4ZrpQqQrRAWKlSMwL//gsBoeV468wrvPt//2V3/HmrUylVIgyaP4hB8wdZHaNI0wJhtXK1CPr7L5iweryT9g5ffvxv1u47ZXUqpVQeZQ33HRsbyy233ELDhg1p0qQJ06dPv6xtYQz3DbB+/Xrmz59fKPnd8ehQGyLSAfgQsAMTjTGX9REXkd7AG4ABNhhjHnQtzwSyRtw6YIyJ9mRWS4VUIGjwfJK/epCRBz/m/c9OEd/zde6PrGR1MqVUHgUFBfHll19St25dDh8+TIsWLbj33nspU6bMhTa5He77WtavX8/mzZvp0KFDoWTPyWNHECJiB8YB9wENgL4i0iBHm7rAy0AbY0xDYFi2t5ONMU1dj+JbHLL4hxI44HvS6nfnH/ZpHP92GBOX77Y6lVIqj+rVq0fdunUBqFSpEhUqVOD48eO5Xn/Xrl3ce++9tGjRgrZt27Jz504Apk2bRqNGjYiMjOTOO+8kOTmZkSNH8vXXX+fr6CM3PHkEcTOw2xizF0BEpgFdgK3Z2jwOjDPGnAYwxsR7ME/R5+OHX6/PyJhfkYFr/sfcRWd46/QoXu7cFLsOGa5UroxeM5rtp7Zfs11Wm9xch7ip3E28ePOLec6yZs0a0tLSqF27dq7XGTx4MBMnTqR27dqsWLGCZ555hoULFzJixAiWLl1KxYoVSUhIIDAwkNdee43NmzczduzYPGfLDU8WiMrAwWyv44CWOdrUAxCRFThPQ71hjMk6oRYgIjFABjDKGDMr5weIyGBgMEC1atUKN71VbDZ8Oo7CUaoSnRb/P8rFPMU/E0bxzoO3EuBrtzqdUiqXjhw5wsMPP8zkyZOx2XJ3siYhIYFVq1bRo0ePC8syMpy3wLdp04b+/fvTq1cvunfv7pHMOXmyQLj7lTdnt2EfoC5wB1AF+E1EGhljEoBqxpjDIlILWCIim4wxl4xPYYz5FPgUnGMxFfYOWMl26xAIrUCrWU9TZs+zPD3+bd4b1J6ywX5WR1OqSMvtb/pZRw6TOkwq9Axnz57l/vvv56233qJVq1a5Xs8YQ1hYmNtrEhMmTGD16tXMnTuXyMhINm7cWJiR3fLkXUxxQNVsr6sAh920mW2MSTfG/AXswFkwMMYcdv3cCywFmnkwa9EU2Qfbg9Op53ucN44PY+i4bzl4KsnqVEqpq0hLS6Nbt24XftvPi7JlyxIREXFhylCHw8GGDRsA2Lt3L61ateLNN9+kbNmyHDp0iNDQUM6dO1fo+5DFkwViLVBXRGqKiB/QB5iTo80s4E4AEQnDecppr4iUFRH/bMvbcOm1i5Kjbjt8HpnLDQGZjE16kVfHTWZjnA4ZrlRRNWPGDJYvX84XX3xx4fbVvNylNG3aND755BMiIyNp2LAhc+fOBeC5556jcePGNG7cmHbt2tGoUSPuuusuNmzYQLNmzTxykdqjw32LSEdgLM7rC58bY/4tIiOBGGPMHBER4H2gA5AJ/NsYM01EWgPjAQfOIjbWGPPZ1T6rSA/3XRhO7iF9clcyzh5jSObz9H1oEHfddH3HhleqqMrPcN+ePMVUVOV1uG+dD8KbnDtG+pfdkePbeCF9MFHRT/Fgy2JycV6pAsjXfBAlkM4HUZyFVsT30Z+R6q0Z4/s/9s15m/fmb784ZPik+y9OgqKUUgWkBcLbBJTC/vD3OBp04xXfqZT5/Q3+Of1P0jIcbDlyhi1HzlidUClVTGiB8EY+/th6fo65+Qke8/mZtlte5bFJK0g0/lYnU0oVI1ogvJXNhtw3Gu5+nS72Pxh88CVeT+zJCUeo1cmUUsWEFghvJgK3PQ9dPqaNfRvv+fyPtxO7kJCUZnUypYq8/Q/3Z//D/a2OUaRpgSgOmj2E9J3GjbY4xvuOYdaqa49Do5QqXNd7uO+ZM2fy7rvvFlp+dzw63Le6juq1J86nGjUy/uL0qik47vgPNh3gT6nrrjCH+87IyMDHx/3XdLdu3Qo/fA56BFGMJEowJyjDPSkL+X33CavjKFUiFXS471tvvZVXX32Vtm3b8tFHHzF79mxatmxJs2bNaN++PfHxzkGvJ06cyLBhzhkS+vXrx9ChQ2ndujW1atW6MFRHQekRBMWoR6UIKbZgGrGPt5Ytom29B61OpNR1d/Ttt0nddu3TrCnbnW1ycx3Cv/5N3PDKK3nOkp/hvsE52N/y5csBOH36NNHR0YgIn3zyCe+//z6jR4++bJ34+HhWrFjBpk2b6N27d6EcYWiBKEYaRpQGRzDph05Q68B3HE7oRqUygVbHUqpEys9w31n69Olz4fmBAwfo3bs3R48eJTU1lXr16rldp2vXrogITZo04dChQwXKnqXEFwiHcVB39gZuPGZ3jgjl7Ww+pNWLpvO2uUxauZ0h95W8QXBVyZbb3/SzjhyqT/my0DPkd7jvLMHBwReeP/3007zyyit07NiRxYsXM2rUZTM3A+Dvf7EfVGENoVTir0EcOL6bZEcqf9QsPreGBt/yKKGSTELMDNIyHFbHUapEKchw3+6cOXOGypUrY4xh8uTJhZAw90p8gaicEULv3xzU25fOjlM7rI5TMIN+cj6qtSKxVG06pS9k4dajVqdSqkQp6HDfOb3xxht069aN22+/nYoVr+8IzjqaK/BbmyaEnUxnyj+b8O9Hp+Echdy7OVb8H7ZF/+Kf4Z/w3tN9rY6jlEflZzRXT55iKqp0NNd8+Gp4JOeDbTSfsZEl+3+xOk6hsDV9kEzxoeGRmew65rkZp5TyVtWnfFmiikN+aIEAxnedQs3nX6bhAfj5qzdJzUy1OlLBBZcno979dLP/zvSVu6xOo5TyQlogXMr3foDMqjfQYV48X238wuo4hcL/5kGUkUTO/TmTxNQMq+Mo5VHF5XS5p+Tnz0cLhIv4+lL95deofAp2f/kJx5Ny3/OxyKp5O6khVenqWMzs2MNWp1HKYwICAjh58qQWiSswxnDy5EkCAgLytF6J7weRXcidd2Br0YSuyzby8e/v8Xr7y3srehWbDb+bB3LLkjf5bMUf9L25arG4AK9UTlWqVCEuLi5PQ1qUNAEBAVSpUiVP62iByEZEqPbyazh69sR36o9sbv4wjcIaWR2rQKRZPxy/vk3UqbmsP3AvLaqXtTqSUoXO19eXmjVrWh2j2NFTTDkENmpIUKeOdFpr+N+Ckd5/yBp6A4467ellX860lbutTqOU8iJaINyo9Pw/sNt8iPxhMz//9bPVcQrM52+DKC9nSN78E6cSi0+PcaWUZ2mBcMO3UiXCBgzkti2Gb+eMIik9yepIBVOnHenBEfSSX5gRc9DqNEopL6EF4grCnngCU6YU9887zhebvXwYcJsd36j+3GbfxOKVa3E4vPy0mVLqutACcQX2kBAihg6j4QH4c+YEjpw/YnWkgmnWDwFuO7+AZbv0Tg+l1LVpgbiKMj17YqtRjb6L0xi75n2r4xRMmWqYWnfxgM8yvlm51+o0SikvoAXiKsTXl0ovvETEKUPmrJ9Zf2y91ZEKxBY1gBs4SeauxcSd9vLrKkopj9MCcQ0hd96B/81R9P4dPlj2bxzGi+dXqHcfmUFh9LEtYeqaA1anUUoVcVogrkFEiHjxJUKTHDSYt43Zu2dbHSn/fPywN3uIu+1/snjNRlIzMq1OpJQqwrRA5EJgw4aUio6mUwxMWTKG82nnrY6Uf80HYMfB3SmLmL9ZJxNSSl2ZRwuEiHQQkR0isltEXrpCm94islVEtojIN9mWDxCRXa7HAE/mzI0Kzw3DbvOhw8KTfLrpU6vj5F/52pjqt9LPbxnfrNxncRilVFHmsQIhInZgHHAf0ADoKyINcrSpC7wMtDHGNASGuZaXA14HWgI3A6+LiKWDCPlGRBA26BFu22JYsfhLDp713g5n0mIglcwx7Ad/Z/vRs1bHUUoVUZ48grgZ2G2M2WuMSQOmAV1ytHkcGGeMOQ1gjIl3Lb8XWGSMOeV6bxHQwYNZc6X8448hZcvQ75cM3lv7rtVx8q9+ZxwBZXjI51e+XqUXq5VS7nmyQFQGsv+aHedall09oJ6IrBCRVSLSIQ/rIiKDRSRGRGKuxzC/9pAQKg4Zwk0HMjm75BdWHVnl8c/0CN8AbJF9uNcew5L1WzmvkwkppdzwZIFwN/FAzjEefIC6wB1AX2CiiJTJ5boYYz41xkQZY6LCw8MLGDd3yvTqhW+tmgxcZuPdlaPIcHjpl2vzAfiYdDpkLmXmn4esTqOUKoI8WSDigKrZXlcBck5rFgfMNsakG2P+AnbgLBi5WdcS4uNDxeHDqXAig5q/7uK7nd9ZHSl/KjbAVPkbA/2X8fXKfd4/rLlSqtB5skCsBeqKSE0R8QP6AHNytJkF3AkgImE4TzntBRYA7UWkrOvidHvXsiIh5I47CGrZkr5/2Ph81f9xJvWM1ZHyRZoPoKojjpD4GGL2n7Y6jlKqiPFYgTDGZADP4Pxi3wbMMMZsEZGRIhLtarYAOCkiW4FfgeHGmJPGmFPAmziLzFpgpGtZkSAiVHhhOEFJmdy1LIH/bfif1ZHyp1F3jF8ID/sv5atV+61Oo5QqYjzaD8IYM88YU88YU9sY82/XsteMMXNcz40x5nljTANjTGNjzLRs635ujKnjehS58bYDGzakdHQ0ndbC4lVT2ZOwx+pIeecXjDTuRUdZzW+bdnPifKrViZRSRYj2pC6A8GFDsdt9eGg5vLv2Xe88j99iAL4mlY6s0MmElFKX0AJRAL4REZQfOIhWm9M4GvM7vx36zepIeVepGdzQhMcCl/P1yv1k6mRCSikXLRAFVP7xx7CXL8fjy335z5rRpGemWx0p75r3p0bGHsqd3crSHfHXbq+UKhG0QBSQPSSE8GefpdZfKVRYt49vtn9z7ZWKmia9MT6BDApcpherlVIXaIEoBGV69sSvVi0eW+7PhPX/42TySasj5U1AaaRhNzqxgjU7D3LwlE4mpJTKR4EQEZuIlPJEGG8lPj5UGP5Pyh5Pps3aRD6K/cjqSHnXYgB+jiQ621fx9Wodn0kplcsCISLfiEgpEQkGtgI7RGS4Z6N5l5A77iCoVSv6/mHn543fsf3Udqsj5U3VlhB2I4ODf2NGzEGdTEgplesjiAbGmLNAV2AeUA142GOpvJCIUPGF4fgnpvHAGl9GrxntXbe9ikCLAdRK3UZ40m5+3qSTCSlV0uW2QPiKiC/OAjHbGJOOm8HzSrqABg0oHR1N+zXp7N+xlkX7F1kdKW+a9MHY/Rgc8jtT9GK1UiVebgvEeGAfEAwsF5HqgM4040b4sKHYbXYeXxnMmHVjSMlIsTpS7gWXR+p35n6znM37j7H1sP4VK1WS5apAGGP+a4ypbIzp6BoeYz+uQfbUpXwjIig3aCBNY88SsDOOL7d+aXWkvGk+gICMs3T2jeGr1XoUoVRJltuL1ENdF6lFRD4TkfXAXR7O5rXKP/Y49vLlGbKiFBM3TiA+yYs6n9W4DcrW4MnQ35n15yHOpXhhxz+lVKHI7SmmR1wXqdsD4cAgYJTHUnk5e0gw4c8+Q8Tu00RuT+PD9R9aHSn3bDZo3p/aSbFUTI/TyYSUKsFyWyCyZnjrCEwyxmzA/axvyqVMz5741a7NE78H8tPO2Ww8vtHqSLnX9CEQO8+W+YMpK/d7191YSqlCk9sCsU5EFuIsEAtEJBRweC6W98vqPBd89AzdNgcxes1oHMZL/shCb4Ab76Nj5q/si09g9V9FZioOpdR1lNsC8SjwEvA3Y0wS4IfzNJO6ipDbbyeoVSu6/57J7rgN/LT3J6sj5V7zAQSknSI6IFbHZ1KqhMrtXUwOnPNC/0tE3gNaG2O86JyJNbI6z/mcS+bxP8sxdt1YktK9ZJyjOndDqSo8WWoF8zcfJf6cF92uq5QqFLm9i2kUMBTnMBtbgSEi8o4ngxUXAQ0aULpLF1qvOI05cozPNn9mdaTcsdmhWT9qn13DDSaeGWt1MiGlSprcnmLqCNzjmgb0c6ADcL/nYhUv4cOGYrPZ+ce6G5i8ZTKHznvJnUHN+iHAP8LW8M3qAzqZkFIlTF5Gcy2T7Xnpwg5SnPnecAPlBg2k1ppD1D7kYEzMGKsj5U6ZqlCnHfdlLObYmUSWbPei/hxKqQLLbYF4B/hTRL4QkcnAOuBtz8UqfrI6zw37oywL9y0g5miM1ZFyp3l/ApKP0TVkm47PpFQJk9uL1FOBVsAPrsctxphpngxW3Dg7zz1LmR2HaX+gDKPXjibT4QVDat94HwRX4KlSv7N853H2nUi0OpFS6jq5aoEQkeZZDyACiAMOApVcy1QelOnZA7/atRmwVNh1fBuzds+yOtK12X2h6YPUOr2CCNtpvlmjkwkpVVJc6wji/as83vNstOInq/Oc7+ETDNpVmf/++V/OpZ2zOta1Ne+PmExejljHjJiDpKR7wZGPUqrArlogjDF3XuWhg/XlQ8jttxN0SyvuWXKa1IRTfLrxU6sjXVv52lDjNtqnLuJMUio/bTxidSKl1HWQ234Q3d087haRCp4OWNw4O8+9gJw9zwvbavPVtq/Yf9YLLv62GEjA+YP0LLtHL1YrVULkZaiNicBDrscE4HlghYjo1KN5FFC/PqW7dKH+L3updM6X99Z6wdm6mzpBYFmeKrWC2IMJbD50xupESikPy22BcAD1jTE9jDE9gAZAKtASeNFT4Yqz8GFDEbudl9ZVZmncUv449IfVka7ONwAi+1Lj+K9U8k3U8ZmUKgFyWyBqGGOOZXsdD9QzxpwCdEaZfMjqPBe2YjttEsL5z9r/kOHIsDrW1TXvjzjSebVyLLNiD3EmWf/qlSrOclsgfhORuSIyQEQGAHNwzk0dDCR4Ll7xVv7Rx7CXL8+TvwWxJ2E3M3bMsDrS1VWoD1Vupl3yfFLSM/lhfZzViZRSHpTbAvE0MAloCjQDJgNPG2MSjTFXnJtaRDqIyA4R2S0iL7l5f6CIHBeRWNfjsWzvZWZbPidvu+UdsjrP+W3ew0PH6zIudhwJKUW83rYYgP+ZPTxY8RBfrdLJhJQqznLbk9oAvwNLgMXAcnONbwYRsQPjgPtwXrPoKyIN3DSdboxp6npMzLY8Odvy6Nzk9EZlevbAr05tuiw4Q3LKOT7e8LHVka6uYTfwL8UTob+z53giK/eetDqRUspDcnuba29gDdAT6A2sFpGe11jtZmC3MWavMSYNmAZ0KUjY4kh8fKg4fDjEHWF4XBNm7JjB7tO7rY51ZX7B0Lgn1Y4upGpgql6sVqoYy+0ppldxziY3wBjTH+eX//+7xjqVcQ7LkSXOtSynHiKyUUS+E5Gq2ZYHiEiMiKwSka7uPkBEBrvaxBw/fjyXu1L0BLdtS9AtrWg2dxdhmYH8Z+1/ivapm+YDkIwU/lVtCwu2HOPYWZ1MSKniKLcFwmaMyT7W88lcrCtuluX81vsR5x1STXCeupqc7b1qxpgo4EFgrIjUvmxjxnxqjIkyxkSFh4dfcyeKqqzOc+bsOf61/SZWHlnJ0oNLrY51ZZWaQkQkdybOI9PhYNoanUxIqeIotwVivogscF1UHgj8BMy7xjpxQPYjgirA4ewNjDEnjTGprpcTgBbZ3jvs+rkXWIrz4nixFVC/PqW7dqXivHW0yKzKezHvkZaZZnWsK2s+AL8TWxlQ/TRT1xwgI9NhdSKlVCHL7UXq4cCnQBMgEvjUGHOtDnJrgboiUlNE/IA+OG+PvUBEIrK9jAa2uZaXFRF/1/MwoA3OqU6LtfChQxC7nWFrwzlw7gBfb/va6khX1rgX+AbxePByjp5NYfE2nUxIqeIm1zPKGWO+N8Y8b4x5zlvu6TkAAByvSURBVBgzMxftM4BngAU4v/hnGGO2iMhIEcm6K2mIiGwRkQ3AEGCga3l9IMa1/FdglDGm2BeIrM5z/r+uoXd6U8ZvHM+J5BNWx3IvoBQ07EbluHnULmX0YrVSxZBc7WKoiJzj8usG4Ly+YIwxpTwVLK+ioqJMTIyXzNJ2FZnnE9nToQOmUgX6dNpDdJ0ujGg9wupY7h1YDZ+3Z0ndf/HIpgYs+cft1AoPsTqVUioPRGSd63rvZa413HeoMaaUm0doUSoOxUlW5znHxq384/ytzNw1k60ni+jBU9WbIfwmbjv3Ez424evVOpmQUsVJrk8xqeunTI/u+NWpTctZOwnzKcPoNaOL5m2vItB8AL5H/+SRuol8G3OQ5DSdTEip4kILRBGU1Xku48BB/nXkb6yPX8+C/QusjuVeZB+w+zEo8DfOpmTw48bD115HKeUVtEAUUcFt2xLc+hYqf7uCSP86jIkZQ0pGEeyQFlQO6kdzw75ZNAz31YvVShUjWiCKKBGhwvDhOM6eZfi2mhxJPMIXW76wOpZ7LQYgKWd4qcZONsadYcPBIj7goFIqV7RAFGFZnef8f1hM9+A2fL75c44mHrU61uVq3AblanFLwlyC/Ox6FKFUMaEFoogLHzYU7Hb6LRMyHZmMXT/W6kiXE4Hm/fE5uJLH6mcyZ8NhEpKKcC9wpVSuaIEo4nwrVqT8I4PIWLSUIX738tPen4iNj3XbdtD8QQyaP+g6J3SJfBBsPvT3X0ZqhoPv1ulkQkp5Oy0QXqDcI49iDwvjtll7CQ8IY/Sa0ThMERv7KLQi1OtA2O7vaVkthK9XH8DhKIK35iqlck0LhBfI6jyX+ucG/pV6D5tPbmbu3rlWx7pci4GQdIJ/VNvNXycS+WOPTiaklDfTAuElyvTojn/dOlT7aimRZRoydt1YktKTrI51qdp3QemqRJ36kXLBfkxZtc/qREqpAtAC4SXEx4cKw4eTfuAgLx6K5HjycSZumnjtFa8nmx2a9cO2dymPN7azeFs8R84kW51KKZVPWiC8SPBttxHc+hb8J8+me8X2TN4ymbhzRexicLN+IMKDvstwGMNUnUxIKa+lBcKLiAgVXniBzLNnGRgTgt1mZ8y6MVbHulTpKlCnHaW3T+euuuWYtuYA6TqZkFJeSQuElwm46SZKd+tG6vSZPB3ek0X7F7H26FqrY12q+QA4d4Qh1fYRfy6VRVuPWZ1IKZUPWiC8UPjQIWC3c+e8Q1QKrsSoNaPIdBShUVTr3QvBFWgSP4vKZQKZslJ7VivljbRAeKGsznOJ8xfyUnB3dp7eyfe7vrc61kV2X2j2ELJrIYObBbBy70l2x5+zOpVSKo+0QHiprM5zNb9cRosKzfnoz4/IcGRYHeui5v3BOOhpW4avXfhqlU4mpJS30QLhpewhwYQPeZbkP//kxcTbSUhN4EjiEatjXVSuFtRsS/CWb7i/UUW+XxdHUloRKmBKqWvSAuHFynR3dp7z/XQ6PWt2IT4pvmjNGdF8ACQc4MnqcZxLzWBOrE4mpJQ30QLhxS52njvAoJ2VsImNA+cOFJ3pSet3hsBy1Iv7gZtuCGXKqv1FJ5tS6pq0QHg5Z+e51iRNmEzXNXA27SwD5w9kx6kdVkcDH3+I7Its/4lHm4ey5fBZYnUyIaW8hhYIL+fsPOecea7DWgc9Y3z568xf9J7bm7dWvUVCisVfyC0GgCOdaJYR7Gdnik4mpJTX0AJRDGR1nit1Jp02O2z82O1H+tzYh+92fkenWZ2Yvn26df0kwm+Eqq3w3/AV3ZtVZu7GI5xO1MmEvFXLST1oOamH1TEKrLjsB3h2X7RAFBPhQ4dgBMqcTKW0f2lebvkyMzrPoF7Zery1+i0emPsA646tsyZc8/5wchePVz9KWoaDb9eVrPGZitOXkSpZtEAUE74VK3KujB/BiRkceeMNHImJ1Ctbj8/af8Z7t7/HmbQzDJw/kBeWv8CxxOs89EXDruBfimr7vuXmGuV0MiGlvIQWiGLkTFk/zpb2JWH6DPZ26UrS2rWICPfWuJfZXWbzRJMn+GX/L3Se1ZmJmyaSlnmdTvX4BUPjXrB1NgObl2b/ySR+233iqqs8MH4lD4xfeX3yKaXc0gJRnIiQEBZA9Slfggj7+w/g2DujcKSkEOQbxDPNnmFW11ncEnELH67/kK6zu7Ls4LLrk63FAMhI4Z7M5YSF+F1zfKZ9fu+xz++965NNKeWWFohiKCgqilqzZlKmzwOcmjyZv7p1J3njRgCqhlblw7s+ZHy78djFzjNLnuGpxU+x78w+z4aKiISIpvjGTuGBqCos2X6MQwk6mZBSRZkWiGLKFhxMxOuvU/WziTiSk9nXpy/xH4zFkeY8rdS6cmt+iP6Bf0b9k/Xx6+k2pxsfrPuAxPREz4VqMQCObaZ/9VMYYOpqHZ9JqaLMowVCRDqIyA4R2S0iL7l5f6CIHBeRWNfjsWzvDRCRXa7HAE/mLC6mPduQac82vGRZSJs21Jozm9JdunBy/Hj29epNyvbtAPjafRnQcABzu83l/pr38/nmz4meGc3cvXM90+O5UU/wDaLirmncfVMFpq09SFqGTiakVFHlsQIhInZgHHAf0ADoKyIN3DSdboxp6npMdK1bDngdaAncDLwuImU9lbW4s5cqRaV33qbKxx+TcfIkf/XqzYlPPsFkOAfPCwsM461b3+Krjl8RHhTOy7+9zID5A9h2clvhBgkoBY26w6bv6d+iPCfOp7Jgy9HC/QylVKHx5BHEzcBuY8xeY0waMA3okst17wUWGWNOGWNOA4uADh7KWWKE3nUntX6cQ6l72nF87Ifs6/sgqXv2XHg/MjySb+7/hhGtR7D/7H4emPsAI1eO5HTK6cIL0XwgpCdya8pyqpYL1J7VShVhniwQlYHsPaLiXMty6iEiG0XkOxGpmpd1RWSwiMSISMzx48cLK3ex5lO2LJXHjKHymPdJP3CAv7p15+SkLzCZzp7WNrHRvW53fuz2Iw/Vf4gfdv1Ap5mdmLp9auHMN1ElCsLrY1s/mX4tq7Pmr1PsPKaTCSlVFHmyQIibZTlPbP8I1DDGNAEWA5PzsC7GmE+NMVHGmKjw8PAChS1pSnXsSK25PxJ8663Ejx7N/gEDSDtw8aJxKb9SvHjzi3zX+Tvql6vP26vf5oG5DxR8/msR58Xqw+vpU+0Mfj42vtKjCKWKJE8WiDigarbXVYBLJgQwxpw0xqS6Xk4AWuR2XVVwPuHhVBn3ERHvvEPq9h3s7dqN01OnXnKBuk7ZOkxoP4Exd4zhXNo5HlnwCC8se4GjiQW4dtDkAbD7U3rbVDo1juCH9YdITNXJhJQqajxZINYCdUWkpoj4AX2AOdkbiEhEtpfRQNZV0QVAexEp67o43d61TF3FpA6TmNRhUp7WERHKdOtKrR/nENS0KUdHjOTgo4+RfuTIJW3uqX4Ps7vO5snIJ1lycAnRs6KZsHECqZmpV9n6FQSVgwbRsHE6/aIqcj41g1mxh/K+HaWUR3msQBhjMoBncH6xbwNmGGO2iMhIEYl2NRsiIltEZAMwBBjoWvcU8CbOIrMWGOlapjzENyKCqp9N5IY3XicpNpa9naNJ+GHmJUcTgT6BPNX0KWZ3nU2bSm3475//peusriw9uDTvt8U2HwApZ2h2bikNIkoxZaVOJqRUUePRfhDGmHnGmHrGmNrGmH+7lr1mjJnjev6yMaahMSbSGHOnMWZ7tnU/N8bUcT3y9muxyhcRoWyfPtSaNRP/m27kyCuvEPfU02TkuAGgckhlPrjzA8bfMx4/ux/PLnmWJ395kr/O/JX7D6txK5Srjaz/kn6tqrP96DnWHyjEu6WUUgWmPanVZfyqVaP65MlUePFFElesYG+nzpz9+efL2rWu1Jrvor9jeNRwNsRvoPuc7oyJGZO73tgizmHAD/xB16rnCfX3ueb4TEqp60sLhHJL7HbKDxpIzZk/4FutGoeee55Dzz9PxulLf8v3tfnSv2F/fuz2I51qdWLSlkl0mtmJH/f8iMNco5d00wfB5kPQpq/p0aIK8zYd5eT5fFzTUEp5hBYIdVX+tWtTY+o3hA8bxtlFi9nbOZpzS5Zc1i4sMIw327zJNx2/ISI4gld+f4X+P/dn68mtV954SAW48T7YMJV+URVJy3QwIybOg3ujlMoLLRDqmsTHh7C/P0HNb2fgU748cU89zeGXXyHz3OUd3BqHN+arjl8xsvVIDp47SJ+5fRixcgSnUq5wj0HzgZB0kjqnf6NVrXJ8s2Y/mTqZkFJFghYIlWsBN91EzW9nUP7vT3Bmzhz2do7m/IoVl7WziY1udbsxt9tc+jXox6xds+g0sxNfb/v68t7Yte+E0lVh3WT6tarOwVPJLN+pveKVKgq0QKg8ET8/KgwbRo2p32ALCuLgo49xZMQIHImXX5gO9Qvlhb+9wHfR39GwfENGrRlF77m9L+2NbbNDs4dh76+0j0glPNRfx2dSqojQAqHyJbBJE2r+8D3lBg4kYdp09nbtRlJMjNu2tcvU5tN7PmXsHWNJSk/ikQWP8M9l/7zYG7tZPxAbfhu/ou/fqvLrjngy00Ou494opdzRAqHyzRYQQMWXXnROcQrsf7g/x0b/B0dKymVtRYS7q9/NrC6zeKrpUyw9uJTOMzszfsN4UkPCoM49EPs1faIqIUBSwo3XeW+UUjlpgVAFdskUp5Mm8Vf3HhemOM0pwCeAJyOfZE7XOdxW5TY+iv2ILrO6sKRmFObcESrF/0a7+hVJPlMP49B/nkpZSf8PVIXikilOk5LY1/dB4seOxbimOM2pUkglxtwxhgntJxBgD2Dorik8WbkKe9d9Sr9W1TGZgaScr3F9d0IpdQkfqwOo4iVritNj74zi5CfjOb90GZVGvUPATTe5bd8qohXfRn/L9O3T+ThmDD3Sd/Ng/EfY/Wpw7lgrunz0OwG+dgL97AT6Oh8B2Z4H+tmd7/vaCfSzOd/P9t6F19me223uRpNXSuWkBUIVuqwpTkPvaceR117nr169CX/6aco/9ijic/k/OV+bL/0a9OO+sg347/c9mLLnB0Jq2jFnG5ASfCPnMyEjDdKTDemZkJEJaRmG9AyDMTYMAggYm/MnAkYuPDdZy40ANnzsNvztdvztPvj7+uBnt+Pv40OAr/Onv48PAa7XAT6+zud+PgT6+BDo60OArw+Bvr4E+foQ6OdDoN/F58F+vgT6+hLo64PdZseGDeOayiTdkQ4GTNZ/5uJPAIdxXPIecLFNjvYO48BhDA6HIdMYDM7XmQ7j/GkycbieOy60B4fD+RmZxkGmMc7XJmsd589Ltm8MmQ7HheWpqX4AfLtpuTNftr/HrOdXGnQx68/h8rcvXZDzfXP5VDAX2lzc5pU+0/12U1IDAMPn6xa6Xa9QebhbT0pKAGLzzIdIcRlBMyoqysRc4S4aZZ2M06c59uabnJ33MwFNmlBp1Dv416p15RUmR7P5zF88EgrJegJUqVxxpNzAlicW5WtdEVlnjIly954eQSiPypriNLRdO46OGMlf3boT/twwyvXvj9jcVIAWA2j03SNEU5HYgECmPLIKk/UbsXFceG7I+i3XzQMHDofz54X2rt+Ks9bNdFzcRvb3L9mG63mmw0FqRgapGZmkZmSQkpFBWkYmqZkZpKZnkJqZSXqm8720zEzSMjNJz/Z8y/HdgFCvXA1EbK4Lf4JNBMl64FouNudynHd+2RBEbIjgWjfrvax2cvE912ubZG3buZ5zG67XODsy2i9534Yt+zYk2zawZdueMPHP6QA80aLvhb+y7CfsJOvVFc7iZb0vcnmDS7Zz2dtycds527saX+nE4aXbcr4Yu/pzAJ5v9egV1ipcnjyp+f6qz/ARz/yirwVCXRelOnYk6G9/48hrrxM/ajTnF/9CxDtv41e16qUNb+oEgeW4O+kcGwOCCPQJtCZwIWo5qQcAsx56w9ogheCr7c6R959q2dHiJAXzvw0fA9C/2V0WJym4j2LHeWzbehCvrhuf8HCqfDyOiHfeIWX7dvZ26crpadMuPX/s4w9NH6RlUhKDZyZbF1YppQVCXV+XTXH6xggOPvb4JVOc0rw/NhtUrpBkXVCllF6kVtYxxpAwfTrH/vMuYrdT8dVXKN2lCyLCqefCKVc6DYLCoFSlSx+hWc8rQ6kI8A+1eleual67BgB0XHyVoc+9RHHZl+KyH1DwfdGL1KpIypriNLh1aw6/8gpHXnqZc4sWEzHiDTbtLkNEWDL1uneGs4fh7CGIWwtJJy/fkH8pN8UjWwEpVRkCy7q78qmUugotEMpyWVOcnvpyCsc/+IC9nTpjbD7sTi5Fvc5jL22cngLnDsPZIxcLx7kjzp9nD0P8djh/FHLOZucTAKGuYpHziCSrsIRUcI4uq5QCtECoIiJritOQtrdx+KWXqbgpgZTzhoNPPoUtMBBbcBASGIgtMAhbUJBzWVAgtqB6SGAktnLBzteBgdj8/RFzDltGApIcj1woIK6icnC1s6hkpuUM4SoiEdmOQCpdWlhCI8DHz5o/JKWuMy0QqkjJmuL0j9aNCUqB9KNHMUlJOJKScCQn40hOhszM3G/QbncWjcBAbEFBSFAQtsDG2IJaYvOzY7M7EB8HNknFRgq2hCRsjvNI5kZs6cuwkYzNx2DzcWDzMYiPwVaqPLZyEUjpyu5PZ4VGgL8OV668nxYIVeSIjw9nSglnSkGzmT9c8p4xBpOWhiMpCZOcfLFwJCXjSEp0LktOxpGYVVCSLrZNvFhkMs+eJSM5ybWec7m5bJjyINfDnROI/QQ2n1jEnukqItkefj5IUCC24BDq2ZJxAMeHdIHLOpJluy6Ss0NXzjY5r6Fc4/XFzmhyhTbXWG67vF0tScEAJ17pjzerKc6/a2/fD3DuS4aHxhfTAqG8iogg/v7Y/P2hbNlC3bbJzMSRnIJJzio6rgKSvchcKCjZlp0/j+PsKRznz+BIPEd6YiImMRnHyTQcaQk40vzBCCcO7CzUvNbwB+D4wbXXaFfUFZf9APAnsEy6R7asBUIVSaly/XtQi92OPSQYQoILdbt5vQ3R7a3nV7odPS/Lsy+7MNqdcbvc4MixjvPnwuibAWg/a5X7z/USC7u2Arx/PyBrX3yp44Fta4FQqohxN07R9b5F90qf5nC9Ywspff3CeEBx2Q+4uC+eoD2plVJKuaVHEKpIGvVgDQC6WRtDqRJNjyCUUkq5pQVCKaWUWx49xSQiHYAPATsw0Rgz6grtegLfAn8zxsSISA1gG7DD1WSVMebvnsyqlKeMeKg+AN49g4JTcdmX4rIf4Nl98ViBEBE7MA64B4gD1orIHGPM1hztQoEhwOocm9hjjGnqqXxKKaWuzpOnmG4Gdhtj9hpj0oBpQBc37d4E/gPk7MaqlFLKQp4sEJWBg9lex7mWXSAizYCqxpi5btavKSJ/isgyEbnN3QeIyGARiRGRmOPHjxdacKWUUp4tEO56b1zosikiNuAD4B9u2h0BqhljmgHPA9+ISKnLNmbMp8aYKGNMVHh4eCHFVkopBZ69SB0HZJ+RvgpwONvrUKARsNTVc/QGYI6IRBtjYoBUAGPMOhHZA9QDdMq4EqJBxGW/DyilrjNPHkGsBeqKSE0R8QP6AHOy3jTGnDHGhBljahhjagCrgGjXXUzhrovciEgtoC6w14NZlVJK5eCxIwhjTIaIPAMswHmb6+fGmC0iMhKIMcbMucrqbYGRIpIBZAJ/N8ac8lRWVfRM6jDJ6ghKlXge7QdhjJkHzMux7LUrtL0j2/Pvge89mU0ppdTVaU9qpZRSbmmBUEop5Za4nZzEC0VFRZmYGL3JSSml8kJE1hljoty9p0cQSiml3NICoZRSyi0tEEoppdzSAqGUUsotLRBKKaXc0gKhlFLKLS0QSiml3NICoZRSyi0tEEoppdwqNj2pReQ4sL8AmwgDThRSHCsVl/0A3ZeiqrjsS3HZDyjYvlQ3xridca3YFIiCEpGYK3U39ybFZT9A96WoKi77Ulz2Azy3L3qKSSmllFtaIJRSSrmlBeKiT60OUEiKy36A7ktRVVz2pbjsB3hoX/QahFJKKbf0CEIppZRbWiCUUkq5pQXCRUTeFJGNIhIrIgtFpJLVmfJLRN4Vke2u/ZkpImWszpRfItJLRLaIiENEvO6WRBHpICI7RGS3iLxkdZ6CEJHPRSReRDZbnaUgRKSqiPwqIttc/7aGWp0pv0QkQETWiMgG176MKNTt6zUIJxEpZYw563o+BGhgjPm7xbHyRUTaA0uMMRkiMhrAGPOixbHyRUTqAw5gPPBPY4zXzCsrInZgJ3APEAesBfoaY7ZaGiyfRKQtcB740hjTyOo8+SUiEUCEMWa9iIQC64Cu3vj3IiICBBtjzouIL/A7MNQYs6owtq9HEC5ZxcElGPDaymmMWWiMyXC9XAVUsTJPQRhjthljdlidI59uBnYbY/YaY9KAaUAXizPlmzFmOXDK6hwFZYw5YoxZ73p+DtgGVLY2Vf4Yp/Oul76uR6F9d2mByEZE/i0iB4GHgNeszlNIHgF+tjpECVUZOJjtdRxe+kVUXIlIDaAZsNraJPknInYRiQXigUXGmELblxJVIERksYhsdvPoAmCMedUYUxX4GnjG2rRXd619cbV5FcjAuT9FVm72xUuJm2Vee2Ra3IhICPA9MCzHGQSvYozJNMY0xXmm4GYRKbTTfz6FtSFvYIxpl8um3wA/Aa97ME6BXGtfRGQA0Am42xTxC015+HvxNnFA1WyvqwCHLcqisnGdr/8e+NoY84PVeQqDMSZBRJYCHYBCuZGgRB1BXI2I1M32MhrYblWWghKRDsCLQLQxJsnqPCXYWqCuiNQUET+gDzDH4kwlnuvC7mfANmPMGKvzFISIhGfdpSgigUA7CvG7S+9ichGR74Ebcd4xsx/4uzHmkLWp8kdEdgP+wEnXolVefEdWN+D/gHAgAYg1xtxrbarcE5GOwFjADnxujPm3xZHyTUSmAnfgHFr6GPC6MeYzS0Plg4jcCvwGbML5/zvAK8aYedalyh8RaQJMxvnvywbMMMaMLLTta4FQSinljp5iUkop5ZYWCKWUUm5pgVBKKeWWFgillFJuaYFQSinllhYIpfJARM5fu9VV1/9ORGq5noeIyHgR2eMaiXO5iLQUET/X8xLVkVUVPVoglLpORKQhYDfG7HUtmohz8Lu6xpiGwEAgzDWw3y/AA5YEVcpFC4RS+SBO77rGjNokIg+4lttE5GPXEcFcEZknIj1dqz0EzHa1qw20BP5ljHEAuEZ9/cnVdparvVKW0UNYpfKnO9AUiMTZs3itiCwH2gA1gMZABZxDSX/uWqcNMNX1vCHOXuGZV9j+ZuBvHkmuVC7pEYRS+XMrMNU1kuYxYBnOL/RbgW+NMQ5jzFHg12zrRADHc7NxV+FIc01oo5QltEAolT/uhvK+2nKAZCDA9XwLECkiV/t/0B9IyUc2pQqFFgil8mc58IBrspZwoC2wBueUjz1c1yIq4hzcLss2oA6AMWYPEAOMcI0uiojUzZoDQ0TKA8eNMenXa4eUykkLhFL5MxPYCGwAlgAvuE4pfY9zHojNOOfRXg2cca3zE5cWjMeAG4DdIrIJmMDF+SLuBLxudFFVvOhorkoVMhEJcU0iXx7nUUUbY8xR13j9v7peX+nidNY2fgBe9uL5uFUxoHcxKVX45romcfED3nQdWWCMSRaR13HOS33gSiu7JheapcVBWU2PIJRSSrml1yCUUkq5pQVCKaWUW1oglFJKuaUFQimllFtaIJRSSrn1/wGnG7p/N35JRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot CV误差曲线\n",
    "test_means = grid.cv_results_[ 'mean_test_score' ]\n",
    "test_stds = grid.cv_results_[ 'std_test_score' ]\n",
    "train_means = grid.cv_results_[ 'mean_train_score' ]\n",
    "train_stds = grid.cv_results_[ 'std_train_score' ]\n",
    "\n",
    "# plot results\n",
    "n_Cs = len(Cs)\n",
    "number_penaltys = len(penaltys)\n",
    "test_scores = np.array(test_means).reshape(n_Cs,number_penaltys)\n",
    "train_scores = np.array(train_means).reshape(n_Cs,number_penaltys)\n",
    "test_stds = np.array(test_stds).reshape(n_Cs,number_penaltys)\n",
    "train_stds = np.array(train_stds).reshape(n_Cs,number_penaltys)\n",
    "\n",
    "x_axis = np.log10(Cs)\n",
    "for i, value in enumerate(penaltys):\n",
    "    #pyplot.plot(log(Cs), test_scores[i], label= 'penalty:'   + str(value))\n",
    "    plt.errorbar(x_axis, -test_scores[:,i], yerr=test_stds[:,i] ,label = penaltys[i] +' Test')\n",
    "    plt.errorbar(x_axis, -train_scores[:,i], yerr=train_stds[:,i] ,label = penaltys[i] +' Train')\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel( 'log(C)' )                                                                                                      \n",
    "plt.ylabel( 'logloss' )\n",
    "plt.savefig('LogisticGridSearchCV_C.png' )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上图给出了L1正则和L2正则下、不同正则参数C对应的模型在训练集上测试集上的logloss。\n",
    "可以看出在训练集上C越大（正则越少）的模型性能越好；\n",
    "但在测试集上当C=100时性能最好（L1正则）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 采用5折交叉验证，用正确率，对Logistic回归模型的正则超参数调优"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, error_score='raise-deprecating',\n",
       "             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,\n",
       "                                          fit_intercept=True,\n",
       "                                          intercept_scaling=1, l1_ratio=None,\n",
       "                                          max_iter=100, multi_class='warn',\n",
       "                                          n_jobs=None, penalty='l2',\n",
       "                                          random_state=None, solver='liblinear',\n",
       "                                          tol=0.0001, verbose=0,\n",
       "                                          warm_start=False),\n",
       "             iid='warn', n_jobs=4,\n",
       "             param_grid={'C': [0.001, 0.01, 0.1, 1, 10, 100, 1000],\n",
       "                         'penalty': ['l1', 'l2']},\n",
       "             pre_dispatch='2*n_jobs', refit=True, return_train_score='true',\n",
       "             scoring='accuracy', verbose=0)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 如果 return_train_score=\"false\"，cv_results_属性将不包括训练分数。\n",
    "grid2= GridSearchCV(lr_penalty, tuned_parameters,cv=5, scoring='accuracy',n_jobs = 4,return_train_score=\"true\")\n",
    "grid2.fit(X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7747395833333334\n",
      "{'C': 0.1, 'penalty': 'l2'}\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print(grid2.best_score_)\n",
    "print(grid2.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gWVfbA8e9JLwRIIAmdUKXX0ARh6YgFFEVdFZZd17K4lt+uorDLqlhA1LUvCorooqgggop0FRelhB567xAInRBIOb8/3glGTCDtzaScz/PM8065d95zs8t7vHdm7oiqYowxxuSWj9sBGGOMKZ4sgRhjjMkTSyDGGGPyxBKIMcaYPLEEYowxJk/83A6gMFWsWFFjYmLcDsMYY4qVFStWHFXVyEv3l6oEEhMTQ1xcnNthGGNMsSIiu7Pab0NYxhhj8sQSiDHGmDyxBGKMMSZPStU1EGNM6ZWSksK+fftITk52O5QiKygoiGrVquHv75+j8pZAjDGlwr59+wgLCyMmJgYRcTucIkdVSUxMZN++fdSqVStHdWwIyxhTKiQnJ1OhQgVLHtkQESpUqJCrHpolEGNMqWHJ4/Jy+/exBGKMMdm47Z2fue2dn90Oo8iyBJID9n8iY0xBKFOmzMX1Pn36UL58ea6//vosyw4dOpQWLVrQqFEjgoODadGiBS1atGDq1Km5+s6VK1cye/bsfMWdHbuIbowxLnjsscdISkrinXfeyfL4W2+9BcCuXbu4/vrrWb16dZ6+Z+XKlcTHx9OnT588x5od64HkwD0HRnLvgRFuh2GMKUG6d+9OWFhYnupu3bqV3r1707p1azp37syWLVsAmDJlCk2aNKF58+Z07dqVc+fO8cwzzzB58uQ89V6uxHogORCo5+kka9m0bB4N2vZ0OxxjTD49/dV6Nhw4dcVyGw56yuRkCLtRlbL864bG+Y4tJ+69914mTJhAnTp1WLx4MQ8++CBz587l6aef5vvvvyc6OpoTJ04QHBzMyJEjiY+P59VXXy3wOCyB5MC70SOpk3AfIbMfJrnJMoJCyly5kjHGeMGJEydYsmQJAwYMuLgvNTUVgI4dOzJo0CBuvfVWbr75Zq/HYgkkB1L8QnmzzEO8cHYkP380jA73veV2SMaYfMhpTyGj5/HpfR28GU6uqCoVK1bM8prI+PHjWbp0KV9//TXNmzdn7dq1Xo3FroHkwKf3deCFxx5mWfj1tD0wmS0rv3c7JGNMKRUeHk7lypWZPn06AOnp6axZswaAHTt20L59e0aNGkV4eDj79+8nLCyM06dPeyUWSyC50GDw6xyVCAK+fpDzyUluh2OMKcauueYabr31VhYsWEC1atWYM2dOjutOmTKFcePG0bx5cxo3bszXX38NwKOPPkrTpk1p2rQpPXr0oEmTJnTr1o01a9bQsmVLu4juprLlK7Cry4s0++FP/PzRk3T482tuh2SMKUbOnDlzcf3HH3/MUZ2YmBji4+N/ta927dpZJpyZM2f+Zl9kZKTXXqRnCSSXmnW9heVrptJm34dsXX0L9Vpc43ZIxhgvKUrXPooiG8LKg/qD3uC4lMNv5oNcOG9TQxtjSidLIHlQLiKS/Z1eoFb6Llb81x4wNMaUTq4kEBGJEJF5IrLV+QzPplwNEZkrIhtFZIOIxDj7RUSeE5EtzrGHCjN+gBY97iCubE9i90xk+9qfCvvrjTHGdW71QJ4AFqhqPWCBs52VD4GxqtoQaAskOPv/AFQHGjjHpng33KzVG/wWJyUMZgwl5cJ5N0IwxhjXuJVA+gGTnPVJQP9LC4hII8BPVecBqOoZVc24d/YB4BlVTXeOJVxavzCUqxDN3qufo07aDuImj3QjBGOMN028zrOYLLmVQKJV9SCA8xmVRZn6wAkR+UJEVonIWBHxdY7VAW4TkTgR+VZE6mX3RSJyr1Mu7siRIwXekJa97mJFWDda7xrPzvVLC/z8xpiSo7Cnc58+fTpjx47Nd9zZ8dptvCIyH6iUxaGcXnX2A64BWgJ7gE/xDF29BwQCyaoaKyI3A+87ZX9DVd8F3gWIjY3VXDQhx2oPepszb7Ul9Yu/kFr/Z/z8A7zxNcaYEqSgpnNPTU3Fzy/rn/KbbrqpYILNhtd6IKraQ1WbZLHMAA6LSGUA5zOrIah9wCpV3aGqqcCXQKtMx6Y569OBZt5qR06ER1ZmV7tnqJe2jeUfP+1mKMaYYiI/07l36tSJESNG0LlzZ958801mzJhBu3btaNmyJb169SIhwfOTOmHCBB555BEA7rrrLh5++GGuvvpqateufXEqlPxw60HCmcBgYLTzOSOLMsuBcBGJVNUjQDcg43HKL53t94EuwBavR3wFra4dwsr1X9B6xzh2b7yZmg1bux2SMSY73z4Bh9ZdudwhZzLCnFwHqdQUrh2dv7hy4dSpUyxatAiA48ePc+ONNyIijBs3jpdffpkxY8b8pk5CQgKLFy9m3bp1DBw4MN89FLeugYwGeorIVqCns42IxIrIBABVTQP+DiwQkXWAAOMz1R/g7H8BuKeQ489Szbv/w1kJJnnaA6Q50ysbY4w33H777RfX9+zZQ69evWjatCmvvPIK69evz7JO//79ERGaNWvG/v378x2DKz0QVU0EumexP45MycC5A+s3w1OqegIocrdGVIiuRlzsSGLjHmPJlFG0v8uGs4wpknLaU8joeQz5xnux5FFoaOjF9aFDhzJ8+HD69u3L/PnzGT066/YFBgZeXFfN/yVhexK9gLXuew+rQjrScutb7NmSt3cYG2NMbpw8eZKqVauiqkyaNOnKFQqIJZACJj4+VB80jmQJIOlzG8oyxmQtP9O5X+qpp57ipptuokuXLkRHRxdglJcnBdGNKS5iY2PVW9MaX2r5jLdps+pJltT/O+1//89C+U5jTPY2btxIw4YNc1epCA9heUtWfycRWaGqsZeWtR6Il8TecD+rg9vTfPPr7NsWf+UKxpiiZ8g3pSp55JYlEC8RHx+q3DWOFPHj1Gf3k56W5nZIxhhToCyBeFFU1Vpsav4kjS6sY/nnL7odjjHGFChLIF7Wpt+DrA2KpcnGf7N/x0a3wzHGmAJjCcTLxMeHqDvfQfHhxJT7bCjLGFNiWAIpBJWq12VD08dpfGENy6e94nY4xpgcGjJ7CENmD3E7jCLLEkghaXPzI6wLbEmT9S9xcPdmt8MxxrggYzr31atX06FDBxo3bkyzZs349NNPf1O2IKZzB1i5ciWzZ88ukPgv5dZkiqWO+PhQ8ffvwPudOfrx/VQatgDxsfxtTGkUEhLChx9+SL169Thw4ACtW7emd+/elC9f/mKZnE7nfiUrV64kPj6ePn36FEjsmdkvWCGqXPMq4hv/jabnV7J8+mtuh2OMcUn9+vWpV8/zHrwqVaoQFRVFbl54t3XrVnr37k3r1q3p3LkzW7Z4JiSfMmUKTZo0oXnz5nTt2pVz587xzDPPMHny5Dz1Xq7EeiCFrM2Av7F+21c0XDuGQ21voFL1um6HZEypM2bZGDYd23TFchllcnIdpEFEA4a1HZbrWJYtW8aFCxeoU6dOjuvce++9TJgwgTp16rB48WIefPBB5s6dy9NPP833339PdHQ0J06cIDg4mJEjRxIfH8+rr76a69iuxHoghczH15fyt7+LL+kkTL4fTU93OyRjjEsOHjzI3XffzcSJE/HJ4ZD2iRMnWLJkCQMGDKBFixYMHTqUAwcOANCxY0cGDRrEhAkTSC+E3xbrgbigau2GLGnwCO03j2HZjLdoe9Nf3Q7JmFIlpz2FjJ7HxD4TCzyGU6dOcd111/Hss8/Svn37HNdTVSpWrJjlNZHx48ezdOlSvv76a5o3b87atWsLMuTfsB6IS9oOHMYG/yY0WPMCCft3uR2OMaYQXbhwgZtuuolBgwZx66235qpueHg4lStXvvhK2vT0dNasWQPAjh07aN++PaNGjSI8PJz9+/cTFhbG6dOnC7wNYAnENT6+voTdNo4AvcCB/95nQ1nGlCKfffYZixYt4oMPPrh4e25u7rKaMmUK48aNo3nz5jRu3Jivv/4agEcffZSmTZvStGlTevToQZMmTejWrRtr1qyhZcuWBX4R3aZzd9mSyc/QfuvLxLUaTeyND7gdjjElVl6mc/fmEFZRZdO5FyNtbhvOJr+G1Fs5iqOH9rgdjjEmk4l9Jpaq5JFbriQQEYkQkXkistX5DM+mXA0RmSsiG0Vkg4jEOPu7i8hKEVktIv8TkWJ7L6yvnx8hA8cRpBfY+6HdlWWMKT7c6oE8ASxQ1XrAAmc7Kx8CY1W1IdAWSHD2/we4U1VbAB8D//ByvF5Vo34LVtX9Cy2TFrPi2/fcDscYY3LErQTSD8h48/skoP+lBUSkEeCnqvMAVPWMqiY5hxUo66yXAw54N1zva3PHSLb41afO8qdJPLzP7XCMMeaK3Eog0ap6EMD5jMqiTH3ghIh8ISKrRGSsiPg6x+4BZonIPuBuYHR2XyQi94pInIjE5WaqgMLm6+dH4IBxhOo5dn/0F7fDMcaYK/JaAhGR+SISn8XSL4en8AOuAf4OtAFqA39wjj0K9FXVasBEINs50lX1XVWNVdXYyMjIPLenMNRs2JoVte+j1ZkfWPmtXbgzxm277x7E7rsHuR1GkeW1BKKqPVS1SRbLDOCwiFQGcD4TsjjFPmCVqu5Q1VTgS6CViEQCzVV1qVPuU+Bqb7WjsLX5/VNs9a1LzNKRHD9y0O1wjDEFqLCnc58+fTpjx44tsPgv5dZUJjOBwXiGngYDM7IosxwIF5FIVT0CdAPigONAORGpr6pbgJ5AiXlXrJ9/AH43v02Zz65l7YdDif3bF26HZIwpYAU5nXtqaip+fln/lN90000FH3wmbl0DGQ30FJGteBLAaAARiRWRCQCqmoZn+GqBiKwDBBjv9Eb+DEwTkTV4roE85kIbvKZW43asrHkPsacXsGruf90OxxhTwPI7nXunTp0YMWIEnTt35s0332TGjBm0a9eOli1b0qtXLxISPIM6EyZM4JFHHgHgrrvu4uGHH+bqq6+mdu3aF6dCyQ9XeiCqmgh0z2J/HJ4L5Bnb84BmWZSbDuS/9UVY67tGsX3MPKr/NIKTrXtSrkK02yEZU2Icev55zm+88nTuyZs8ZXJyHSSwYQMqDR+e61jyMp07eCZjXLRoEQDHjx/nxhtvREQYN24cL7/8MmPGjPlNnYSEBBYvXsy6desYOHBgvnso9iR6EeUfEAj93qKcnmbLhw+6HY4xxgvyMp17httvv/3i+p49e+jVqxdNmzbllVdeYf369VnW6d+/PyJCs2bN2L9/f75iB5vOvUir0+xqliz7A+33vcfqBVNo0f32K1e6gtve+RmAT+/rkO9zGVNc5bSnkNHzqPnRhwUeQ16nc88QGhp6cX3o0KEMHz6cvn37Mn/+fEaPzvrJhsDAwIvrBTEPovVAirhWdz/PTp8Yqvz4JCePH3U7HGNMAcjPdO5ZOXnyJFWrVkVVmTRp0pUrFBBLIEVcQGAQqTe+SYSeYPMke/GUMSVBfqdzv9RTTz3FTTfdRJcuXYiOLrzrpTadezGxZPxDtN8/ibVd3qNZ11vyfB4bwjKlVV6mc/fmEFZRZdO5l0At7x7NLp/qRP/wOKdOJLodjjGlQs2PPixVySO3LIEUE4FBIVy4/k0q6jE2ffiw2+EYY4wlkOKkfqvfsbzy72l77CvWLcrq4X1jzOWUpiH7vMjt38cSSDHTYtCL7JUqVFz4N86cOu52OMYUG0FBQSQmJloSyYaqkpiYSFBQUI7r2HMgxUxQSBnOXvs69b+5leWTHqHdXwvvlj1jirNq1aqxb9++XE0ZUtoEBQVRrVq1HJe3BFIMNWjbkyUrbqP94SnEL/6KJh1vcDskY4o8f39/atWq5XYYJYoNYRVTzQe9xD6pTMT8v3H29Em3wzHGlEKWQIqp4NAwTvV+lUrpCcR/+H9uh2OMKYUsgRRjjdr3YXnUANodmcqGJbPdDscYU8pYAinmmg5+hQMSTdk5j3Du7Gm3wzHGlCKWQIq5kDLlON7jFarpQdZ8+He3wzHGlCKWQEqAxh2vZ2mF/rQ99Cmbls1zOxxjTClhCaSEaDz4VQ5LRUJmP0xy0hm3wzHGlAKuJBARiRCReSKy1fkMz6JMVxFZnWlJFpH+zrFaIrLUqf+piAQUfiuKljJlw0ns9hI10vez6qNhbodjjCkF3OqBPAEsUNV6wAJn+1dU9TtVbaGqLYBuQBIw1zk8Bvi3U/848KfCCbtoa9K5P8sibqDtgclsjlvodjjGmBLOrQTSD8iYg2MS0P8K5W8BvlXVJBERPAllai7qlxoNBr3GUYkgcNZDJJ9LcjscY0wJ5lYCiVbVgwDOZ9QVyt8OfOKsVwBOqGqqs70PqJpdRRG5V0TiRCSuNMyBU7Z8BQ53eZGY9L2s+uhJt8MxxpRgXksgIjJfROKzWPrl8jyVgabAnIxdWRTLdnpNVX1XVWNVNTYyMjI3X11sNet6C8vK96XN/g/ZuvpHt8MxxpRQXptMUVV7ZHdMRA6LSGVVPegkiITLnGogMF1VU5zto0B5EfFzeiHVgAMFFngJcdXgNzj2Wht8Zz7IhYZLCQj0TNE8MvExp8T/3AvOGFMiuDWENRMY7KwPBi73dqQ7+GX4CvVM5v8dnusiOalfKpULr8jBa16gdvouVnw0wu1wjDElkFsJZDTQU0S2Aj2dbUQkVkQmZBQSkRigOvDDJfWHAf8nItvwXBN5rxBiLnaad7+d5eV6Ebt3ItvX/uR2OMaYEsaV94GoaiLQPYv9ccA9mbZ3kcUFclXdAbT1YoglRv1Bb3LyjbYwYygpDZa4HY4xpgSxJ9FLuHIVotl79XPUSdtB3OSRbodjjClBLIGUAi173cWKsG603jWexNRgt8MxxpQQlkBKidqD3uaMhBKdeoC0dLejMcaUBJZASonwyMrsavcMV/ns41CK9UKMMflnCaQUaXXtEOLS69ORNezYvsXtcIwxxZwlkNLGPxg/0tj72WOk2liWMSYfLIGUMsG+6ayTenQ5/z3ffD3N7XCMMcWYJZBSKNjPh2N+UdRfOYpth066HY4xppjKdQIRER8RKeuNYEwh8fHBt89zNJTdLPzvaNLSs52L0hhjspWjBCIiH4tIWREJBTYAm0XksSvVM0VXuda3cqRiWwaensTH3610OxxjTDGU0x5II1U9hefFTbOAGsDdXovKeJ8IFW99jTA5h98Pz7Pz6Fm3IzLGFDM5TSD+IuKPJ4HMcKZWt3GPYuiZCmN5psJYACS6Eckt/8RtPgv4zydfkG5DWcaYXMhpAnkH2AWEAotEpCZwyltBmcIT2usfXAgI59Yjr/PhTzvdDscYU4zkKIGo6uuqWlVV+6rHbqCrl2MzhSG4PIF9nqaNzxY2zH2PPYn2HnVjTM7k9CL6w85FdBGR90RkJdDNy7GZQiIt7uJCdAv+LpMZ+fnPNpRljMmRnA5h/dG5iN4LiASG4LwEypQAPj4E3PAyUXKc9vveZ/KyPW5HZIwpBnKaQMT57AtMVNU1mfaZkqBaLNriTu7xm82UWQvYe8yGsowxl5fTBLJCRObiSSBzRCQMsImUShjp8RQ+AcE8KR/w5LS1eF4/b4wxWctpAvkT8ATQRlWTgAA8w1h5IiIRIjJPRLY6n+FZlOkqIqszLcki0t85NllENotIvIi879xibPKrTBQ+XYfTSdYQvHMOU5bvdTsiY0wRltO7sNKBasA/ROQl4GpVXZuP730CWKCq9YAFzval3/mdqrZQ1RZ4LtgnAXOdw5OBBkBTIJhM71E3+dT2z2hkA54N/pix36xh/4lzbkdkjCmicnoX1mjgYTzTmGwAHhKRF/Lxvf2ASc76JDwPKF7OLcC3Tu8HVZ3l3E6swDI8yc0UBF9/5NoXiU47xGCdyZNfrLOhLGNMlnI6hNUX6Kmq76vq+0Af4Lp8fG+0qh4EcD6jrlD+duCTS3c6Q1d3A7PzEYu5VO0u0KgfQ/1msm3LRj5fsc/tiIwxRVBuZuMtn2m93JUKi8h85xrFpUu/3AQoIpXxDFXNyeLw28AiVf3xMvXvFZE4EYk7cuRIbr66dOv1LL4+wsvlP2fU1xs4dDLZ7YiMMUVMThPIC8AqEflARCYBK4DnL1dBVXuoapMslhnAYScxZCSIhMucaiAw3Zl/6yIR+ReeZ1L+7wpxvKuqsaoaGxkZecWGGkf5GkinR+mQ/COt0tYyfLoNZRljfi2nF9E/AdoDXzhLB1Wdko/vnQkMdtYHAzMuU/YOLhm+EpF7gN7AHc4FfuMNHR+C8jV4tewnLNp0gOmr9rsdkTGmCLlsAhGRVhkLUBnYB+wFqjj78mo00FNEtgI9nW1EJFZEJmT6/higOvDDJfXHAdHAz84tviPzEYvJjn8w9H6B8LPbGR65mKe/2kDCKRvKMsZ4+F3h+MuXOabkcT4sVU0EumexP45Mt+Sq6i6gahblrhS3KSgNroM63Ri89xPGp7RixJfxvHt3a0RsIgJjSrvL/hCrqs24CwyZ7XlmcmKfiS5H4gIR6DMG3/90YGL1WfTZcBsz1xygX4vf5HVjTCmTo/+SF5Gbs9h9Elinqpe7AG5Kgsj60P4BGvz0BrdW+h1PzVzP1XUqEhkW6HZkxhgX5WYqkwnAnc4yHs/dT4tFxF5tWxp0fhzKRDMqYBJJF1L418x4tyMyxrgspwkkHWioqgNUdQDQCDgPtAOGeSs4U4QElYWezxCUsJr/NN7MrHWH+GbtQbejMsa4KKcJJEZVD2faTgDqq+oxICWbOqakaXYbVGtL131v076KHyNnxJN45rzbURljXJLTBPKjiHwtIoNFZDCe5zgWiUgocMJ74ZkiRQT6jkXOHuXtqnM5lZzCU19tcDsqY4xLcppAhgITgRZASzwTIA5V1bN2p1YpU6UFtB5MRPxEnmrnw1drDjA7/pDbURljXJCju7BUVUXkf8AFPM9/LFOb16JY+vS+Dvk/SbeRsP5L7jj2Fh9Xfox/fBlPu1oRhIcG5P/cxphiI6fTuQ/EM236LXjmploqIrd4MzBThIVWgG7/wGfXIsa13s+JpAs887UNZRlT2uR0CGsEnrcRDlbVQUBb4J/eC8t4y5DZQy4+GJkvrYdAdBOqL3+OhztXY/qq/czfcPjK9YwxJUZOE4jPJQ8MJuairimJfP3g2hfh5F7+4v8VDSqFMXz6Ok4m2U15xpQWOU0Cs0Vkjoj8QUT+AHwDzPJeWKZYiOkITW7B96fXeK13OIlnLzDqGxvKMqa0yOl07o8B7wLNgObAu6pqDxAa6DUKfPy4as1oHuhSh6kr9vHdZpvdxpjSIMfDUKo6TVX/T1UfVdXp3gzKFCNlq0Dnv8Omr3koZg/1o8vw5LR1nEq2oSxjSrorvQ/ktIicymI5LSKnCitIU8R1GAoRtQmY+yQv3dSQhNPJPP/NRrejMsZ42WUTiKqGqWrZLJYwVS1bWEGaIs4vEPqMhsStNNs/hXs712HK8r0s2mLvoDemJLM7qUzBqN8b6vWGH8bwSLsw6kSG8uQX6zhzPtXtyIwxXmIJxBScPi9A2gWCvn+GF29pzoGT53hhlg1lGVNSWQIxBadCHejwIKydQmvZwp861mLy0j38tO2o25EZY7zAlQQiIhEiMk9Etjqf4VmU6SoiqzMtySLS/5Iyb4jImcKL3FzRNX+DsCrw7WP8rUddalUM5fFpazlrQ1nGlDhu9UCeABaoaj1ggbP9K6r6naq2UNUWQDcgCZibcVxEYoHyhRHs7W+s5/Y31hfGVxV/gWU8z4YcXENw/GRevKUZ+0+c48XZm9yOzBhTwNxKIP3wTAmP89n/MmXBM4njt6qaBCAivsBY4HGvRWjyrskAqNkJFjxDmygY3CGGST/vZsmORLcjM8YUILcSSLSqHgRwPqOuUP524JNM2w8CMzPOcTkicq+IxIlI3JEjdltpoRCBa8dA8gn47nke73MVNSJCGDZtLecupLkdnTGmgHgtgYjIfBGJz2Lpl8vzVAaaAnOc7SrArcAbOamvqu+qaqyqxkZGRua2GSavKjWB2D9B3HuEHNvImAHN2J2YxNg5m92OzBhTQLyWQFS1h6o2yWKZARx2EkNGgrjc5EkDgemqmjE3RkugLrBNRHYBISKyzVvtMPnQdTgElYdZj9OhdgR3t6/JxJ92ErfrmNuRGWMKgFtDWDOBwc76YGDGZcreQabhK1X9RlUrqWqMqsYASapa12uRmrwLiYDuI2HPTxA/jSeubUDV8sE8PnUtySk2lGVMcedWAhkN9BSRrUBPZxsRiRWRCRmFRCQGqA784EKMpiC0GgSVm8PcfxBKMmMGNGPH0bO8Mm+L25EZY/LJlQSiqomq2l1V6zmfx5z9cap6T6Zyu1S1qqqmX+ZcZQojZpNHPr7Q9yU4fRB+fImOdStyR9saTPhxByv3HHc7OmNMPtiT6Mb7qreF5nfAz29B4naG921ApbJBPPb5GhvKMqYYswRiCkePp8E3EGY/SViQPy8MaMb2I2d5bcFWtyMzxuSRJZAcCD2VQpmTF0hPTnY7lOIrLBp+Nwy2zoEtc+hSP5KBsdV4d9EO1u474XZ0xpg8sASSA2d8LhBx9Dybu/6Oo+PGkXbypNshFU9t74OK9WH2E5B6nhHXNSKyTCCPfb6W86k2lGVMcWMJJAemdfbjqd/7sjL8FEdefY2Nv+vMvuefJeXQIbdDK178Ajwvnjq2A35+k3LB/rxwc1M2Hz7NmwvtUR5jihtLIDlwy8pAbl0dTMqLj/Pvh2vwU+0UTnw0mc3duxP/f/eTvM1+/HKsbndocD0seglO7qdrgyhublWVt7/fTvz+nPfsbnvnZ25752cvBmqMuRJLIDkUdl74Q5M/8M79s2n7n4/5ZnRfFrbyJXXuD+y8/gZ+uvtGDi21x1VypPdzoOkw758AjLy+ERGhATw2dS0XUrO9Y9sYU8RYAsklEaFFVAuG9XuZP73/M7s/GM6iXpXwW7eV44PvZ8ENV7Nk2lukptn7L7IVHgMdH4b4abDrf5QPCeD5m5qy8eAp3v7eenPGFBeWQPIh1D+UfrF3c9/r3xEx6/NBUxYAABxSSURBVHPW39mOoMMnKTfiTb7v1pLPXn+QPcd3uR1m0dTxEShXHWY9Dmmp9GwUTb8WVXhz4TY2HjzldnTGmBywBFJA6lZpwi3//IA2i5Zz7O93ESQBNH17ATv7XMubI/oya8N0klPtNuCLAkI8Q1kJ62HFRACeuqEx5UP8eWzqGlLSSs9Q1pDZQxgye4jbYRSIdhMH0G7iALfDKBAlpS3ebIclkAIWEBRCx3tG0HHhckL//RyBlarQfdpOKtw5nH8/2IGx80eyMXGj22EWDQ1vhFpdYOGzcDaR8NAAnu3fhPj9p3jnh+1uR2eMuQJLIF4iPj7UuPZm2s9YQPX/fkhwi+b0/z6ZHo9+zjcPD+DPH/bnk02fcPJ8KX6mRASufREunIGFzwDQp0llrmtWmdcXbGPL4dMuB2iMuRxLIIWgTGwbmk2cQq2ZMwjv05e+q4SHR2/m5D9Gcfdbv2PYomEsPbiU9OznjCy5ohp4HjBcMQkOrALgmRsbUybIj8c+X0NqKRrKMqa4sQRSiILq1ydm7MvUmzePinfdTZdtAYx5N5kWL33L2Pf+RN8v+vLOmnc4dLaUPaD4u2EQWtFzQT09nQplAnmmX2PW7DvJ+B93uh2d1204eIoNduOAKYYsgeTAlL82ZspfGxfY+fyrVKHS8OHU++47Kv71QVolhDLqv2n8fcIxfv78dfpM7cUD8x9g3u55pKSlXPmExV1QOejxFOxbBms/BeC6ppXp07gS/56/hW0JZ1wNzxiTNUsgLvILDydy6FDqLVxA9IgRxJwPY9jUdMZ/FEb5hat5bMGj9Jjag7HLx7L9RAm/qNz891A1FuaNhORTiAij+jchJMCXx6auIS1d3Y7QGHMJSyBFgE9ICBF330WdObOpMvZFKoRUZNAXJ/jvxLIMjo9g2trJ9J/Rnztn3cm0LdM4m3LW7ZALno8P9H0Rzh6BH8YAEBkWyNM3NmbVnhNMXPzroaxdAS+xK+AlNyI1xjgsgRQh4u9PuRtuoNaML6n+zjjCatahw+eb+ODdQP69rR16/CRP/fwUXT/ryj8X/5NVCatQLUH/ZV61NbS8C5aOgyObAbixeRV6NIxm7JzN7DhiQ1nGFCWuJBARiRCReSKy1fkMz6JMVxFZnWlJFpH+zjERkedEZIuIbBSRhwq/Fd4jIpTp0oWa//2Imp98TGhsG6p+vpiRY/fzyebu3BLWhbm75jLo20H0m9GPifETOXruqNthF4zu/wL/UM+U76qICM/f1IRAPx+GTVtLug1lGVNkuNUDeQJYoKr1gAXO9q+o6neq2kJVWwDdgCRgrnP4D0B1oIGqNgSmFErULghp2ZLqb71J7W++pmzfvvh+tZAbnvyWz1d2ZEzlBygXUI5XVrxCz8978vDCh/lh7w+kpmc/D9ftb6zn9jfWF2ILcqlMJHQdDtsXwqZvAIgqG8S/bmjM8l3HmfTzLlfDM8b8wq0E0g+Y5KxPAvpfofwtwLeqmuRsPwA8o+p5cEJVE7wSZRESWKcOVZ5/jrrz5hJx992c++4Haj30Bs99Gcz0Gs9xV8M7WX1kNQ8ufJDeU3vz2srX2HNqj9th502beyCqEcx5ElLOAXBzq6p0vSqSMbM3sTuxBF4DMqYYciuBRKvqQQDnM+oK5W8HPsm0XQe4TUTiRORbEamXXUURudcpF3fkyJF8B+42/0qViH5iGHW/W0jkIw+THL+elAeGMeClOL4MH86rnV+hQYUGvB//PtdNv44hs4fw1favOJd6zu3Qc87XD64dAyf2wOLXAc+w3vM3N8Xfx4fHp66lJF36Maa48loCEZH5IhKfxdIvl+epDDQF5mTaHQgkq2osMB54P7v6qvquqsaqamxkZGRemlIk+ZYrR8X776fuwgVEj/wnqYmJHHroEWKGvspzJ7oz58ZveKjlQxxOOszw/w2n+2fdeXbJs+wLT0cpBr++tTpDo/7wv1c8iQSoXC6Yf1zfkKU7j3HuRAOXAzTGeC2BqGoPVW2SxTIDOOwkhowEcbkhqIHAdFXN/ETdPmCasz4daOaNNhQHPkFBRPz+99SZ/S1VXn4JCQjg4IgRnO53FzfH+TGz16e83/t9ulTvwpfbvuTVnud5qfd5xq8dz77T+9wO//J6PQsIzBlxcdfA2OpcU68iZ460Ie1CGfdiM8a4NoQ1ExjsrA8GZlym7B38evgK4Es8F9YBugBbCjS6Ykj8/Ch33XXUmv4F1cePJyAmhoQXX2R7tx7UnLKYUY3+zsKBC7l5hT8hF4TXV73OtV9cy12z7uLjjR+TeC7R7Sb8VvnqcM3fYONM2PE94BnKGj2gGYhy6nAnzpy3F3cZ4xa3EshooKeIbAV6OtuISKyITMgoJCIxeO62uvRdsaOBASKyDngBuKcQYi4WRIQy13Si5qQPiPnsU0LbtSPxnXfZ1q07Z1/4N9ds8mHod4HMGTCHR1o9wrnUc7yw7AW6f96d++fdz8ztMzlzoQg9b3H1Xz1vMJz1ODjTulQtH0yZyOVcSKpC86fn0u/N//H8rI0s2HiYk+dKwdQvxhQRfm58qaomAt2z2B9HpmSgqruAqlmUOwFc58UQS4TgZs2o9sbrnN+xk2MT3+fk1GlUSUnhfJAvfhM+47Y2bRjS4w62n9/Ptzu/ZdbOWYz43wgCfQPpUq0LfWv35Zqq1xDgG+BeI/yDoPcLMOUOWPYudBgKQMPQWZyqvIGeDcewdMcxPli8i3cX7UAEGlUuS7taFWhXO4J2tSIoH+Ji/MaUYK4kEFO4AmvXovKoUVR88K+s6N+ToKRUEsdPIHHcO+DnR3Djxtzetg1DYp9keyt/ZiV8z5xdc5i7ey5h/mH0jOlJ31p9iY2OxdfHt/AbcNW1ULcHfD8amt4KZTw37ZUN3sXfel0FQHJKGqv2nGDJjkSW7kxk8tLdvO9Mf9KgUhjtakXQrnYF2taKoGKZwMJvgzElkCWQUsQ/OoqTEYGcjAikx6ffcW7VKpKWLydp+XISP5gE4ycQ6OPDHY0aMST2WvbWCWNWuV18u3M2X2z9gsjgSPrU6sN1ta6jUYVGiEjhBC4CfcbA2+1h/tPQ/63fFAny96VDnQp0qFMBgPOpaazZe5KlOxJZuvMYn8XtY9LPuwGoG1XmYkJpXyuCqLJBhdMOY0oYSyCllG+ZUMpc04ky13QCIP3cOc6tXu1JKMuWc3Lyx4SlpHCbCHdfVZ+jDaL5Keok049/wkcbPqJm2Zr0rdWXvrX6ElMuxvsBV6wLHf4Ci1+D2Cu/PzzQz5e2tSJoWyuCvwIpaems3XeSZTuPsXRnIjNWH2DyUs/twbUqhjoJJYJ2tSpQpXywlxtjTMlgCcQA4BMcTGiHDoR26ABAenIy59asvdhDKfftUvqcP08fILlmNBtrnOX7Cm8zucbbVKvRmL61+tInpg/RodHeC7LzY7DmU5j1GH/47DwgcG/Oqvr7+tC6Zjita4bzwO/qkJqWzoaDp1i6w5NQZq07yJTlewGoHhHsuYZSK4L2tStQLTy48HpbxhQjlkBMlnyCgght15bQdm0BSL9wgeT4eJKWeRJKqxWraJmUBsCRqM2srrqOf1UfS0BsS7q06E+Pmj0oF1iuYIMKDIOez8D0e6kWXZ59h0PzfCo/Xx+aVStPs2rl+XPn2qSlK5sOnWLJjmMs3ZHI/I2HmbrC85xMlXJBtKtd4eKwV0yFEEsoxmAJxOSQT0AAIa1aEdKqFdx/H5qSQvKGDSQtX07o8uVExsXRc1USzIzjYHgck2v8i/QWDbiq+wA6tepPsF8BDQs1Gwhx73PVhWUcOlpwQ02+PkLjKuVoXKUcf+pUi/R0ZUvC6Ys9lB+3HmH6qv0ARJcNpO3FHkoEdSLLWEIxpZIlEJMn4u9PcPPmBDdvToV77kFTU0netJmk5cvQxQupsHItAWvWw6T1LC8/ilMNqxPdqRtNegwkqEZM3n9wRaDviwSM60y9mt57j7iPj9CgUlkaVCrL4KtjUFW2Hznj6aHs9PRSvlpzAICKZQJoWyvi4q3D9aPC8PGxhGJKPksgpkCInx/BTRoT3KQxFYYMQdPSSNq0iS0Lp3P+5++JWr2XMj9/wK6xH5AUEUJg65ZU7dSLkLZtCIjJZUKp3Jw9h0KJqXwW3usN4TU9DxtmLOVrQlhlz1sOC6p9ItSNCqNuVBh3ta+JqrIrMeniXV5LdyQya90hAMJD/GkTE3Fx2Kth5bL4Xiah/GvyRs/Kle8NMKZIsQRSykz5a2MAenv5e8TXl9DGjWnZuDH89R+cT0lm6U9T2fr9l/is2USDnxZzaN5iT+GKEZRt246QNm0IadOGgDp1rphQNu8qiyrExPjBrsWw9jPIPEmkbyCUr+EklSwSTFDZ/LVPhFoVQ6lVMZTb29YAYO+xJOc5FM+w19wNhwEoG+TnJJQI2sZEUKOikHj+CAlJCSQkJfBjU+FcAMQvfzFfMRWU9HQlTZXUNCUtXUlN93ympaeTmn7p/nRS0yEtPZ20dCVsfwqKMODT37zip9gJ238BoNi3JaMdm48c4KrIKgV6bilRr0S9gtjYWI2Li3M7DFcNme35z9yJfSa6FsPZlLMs3L2AxUunciFuFQ12p9Fsvy/lTnrmtfKNiCAkNtaTUNq2IbBePeSS3sSsHo0A6Dt/g2dH6nk4uQ+O74Tju+D4bufTWT9/8tdBBEdkSiqXJJiy1TxTyudSSnoKiecSOZx0mISkBLYe3c+6w7vZcfwgR5ISOK/HEf+TiM9vp1sJSFH8Q7KZHFJ/SY3qbGT+V5vxTzhjluXM/6Qve+ziufVX2/mnnqHG4i7jj1Xc2+K04z9dP+KaWo3ydAoRWeHMfv4r1gMxhS7UP5Qb6t7IDXVvJPHmRObunstb27/h4NbVNNqjdEoQrlq1lNNzPS+g9C1XjuDYWELaeJJKUIMspnL3C4QKdTxLVs4dz5RQdv2SYA6s8kzWmPktjuIL5apdTChavgZnylUmISiMwwGBJKSeI+Gcp/eQkSwSkhJIPJf4m6ny/X38iQqJoklEFOX8a5JyoSzHTwWz/6g/BxMDSU8ti//5IMLPHSOkTiOSU9JITkknOTWN5JQ0zqem5/ndJ34+QpC/L0H+PgT6eT492866ny9BAb4EZuz/VZlf9gVmrud3yTkylQn08+Hbnp4e7sXEXoz95j9SiqmMdlzzx7wlj8uxBGJcVSG4Anc0uIM7GtzBvtP7mL1rNp/s+IZtJ7ZR6WQAN5yqTbvDYcjGTZxZsAAAn7AwolOU5AA4Nnky4ueP+Pkifn7g64f4+SF+vuDnh/j6If5+iK+z7VcXqdAAonwRP3/SfODEheMcPbmdo6d2kHB6DwnnDnD4/DEOp+7lcMI2DicK57K4nlIOX6L8QogKjKBBWB2iqnUjKrwO0WHViAqJIiokivDA8GyH446eOc+yncf49JXXORIUSdWIkEt+pH/5sQ689Ifb75If8YtJ4pd9/r5uzZVqSgtLIKbIqBZWjXua3sM9Te9hy/EtzNoxi693zmJ89W0EtQviujI96HOyOtW3nSZ5+lQikuHwqGcL5Lv9gcrOcql0X0F9FfH1QXwEHwEfSceHC0ASIgmIbEJ8FATEPwAJCORMYAhnAkOQoFAkqAwEl/Ws+/mDk/Ba+foRse5LFKhT3wfx9UX8s0mEWSRGfH0vSaBOYvTzJT0jof4qgTqLry/4+3u+z9eF+c1MiWAJxBRJ9cPrU791fR5q9RBrjqzhmx3fMHfXXKbxPWUblaVBug/Nd6XT9PEXOHr6EIlnjnD87BGOnz3K8bNHOXH2GCkpyfim4yyKbxqE+YYQ4VeOcP9yhPuGUc6vDOV8QynrW4ayPiGU8Q0hGH9IS0dTUyAtDU1JRdPSPNupaWhqKpqWCqmpnvXks+i505B82rN+PgnOnyP91HE09TCaLp5rDemg6oOKH+DrWVeh7DnPRY3E8eMhPb3w/9gi4OvjJBZPEvMkIt9fEkym42Ta70lMzrqfL1VOpoHA/gf+iKfRysULKxfXNdMFmMzH0zOVy6Z85jq/2ccv+zPOlV0MWZ0v47hzrvop5wHYP/CaAvgjuyejHSnb1+Ffp2mBntsSiCnSfMSHllEtaRnVkmFth7HkwBJm7ZzF3LNfsay+H6z+JwB+4kdkSCRRFaOIqtGEuiHRF4eRokKiiA6JJjIksuAeaMyplGQ4ufeS6y+7frkGc+H0r4pf/H1LB1UB51PT+SUROduecuL5zUu/pGxGXec4F8sBmcpcdvvSY+czH/PEkp75+xXC8JRJjvuxsP7COZRpGPFXQ4rZ79dUzxBgcnIRfNlaLmS0Q88W/HNTlkBMseHv48811a7hmmrX0O7ZGeyJhB6vTSUqJIqIoAh8pAiO+fsHQcV6nuVSqpB0jMV3tiPQP43YZ/9z8eesuN73E/fPB1AV2jz/jucHWXw8NyWIzy+LT8a2XOZYTo77XOFY/v6KJe0ier1mHQv83JZATLEUkAZ1D0GjCgV/Z0mhEYHQCpw847zw6qo+7sZTABKOOT28ej3dDcQUCksgplg6LzblujFuK4J9fmOMMcWBKwlERCJEZJ6IbHU+w7Mo01VEVmdakkWkv3Osu4isdPb/T0TqFn4rjJtCAnwJCSgZt59+MDCADwbae9tN8ePWENYTwAJVHS0iTzjbwzIXUNXvgBbgSTjANmCuc/g/QD9V3SgifwH+AfyhkGI3xmTj6TsbAtDX5TgKQklpizfb4VYC6Qf8zlmfBHzPJQnkErcA36pqkrOtQMZseOWAAwUfoinKCmtSSGNM9txKINGqehBAVQ+KSNQVyt8OvJJp+x5gloicA04B7bOrKCL34rz4tEaNGvkK2hhv2O2fzfxdxhRxXrsGIiLzRSQ+i6VfLs9TGWgKzMm0+1Ggr6pWAyby6+TyK6r6rqrGqmpsZGRkXppijDEmC17rgahqj+yOichhEans9D4qAwmXOdVAYLqqpjh1I4HmqrrUOf4pMLug4jbGGJMzbt3GOxMY7KwPBmZcpuwdwCeZto8D5USkvrPdE9hY4BEaY4y5LLeugYwGPhORPwF7gFsBRCQWuF9V73G2Y4DqwA8ZFVU1VUT+DEwTkXQ8CeWPhRq9McYYdxKIqiYC3bPYH4fnAnnG9i6gahblpgPTvRiiMcaYK7An0Y0xxuSJJRBjjDF5YgnEGGNMnlgCMcYYkyeWQIwxxuSJJRBjjDF5YgnEGGNMnlgCMcYYkyf2SltjXLZ0yDS3Qygw1paix5vtsB6IMcaYPLEEYowxJk8sgRhjjMkTUVW3Yyg0sbGxGhcX53YYxhhTrIjIClWNvXS/9UCMMcbkiSUQY4wxeWIJxBhjTJ5YAjHGGJMnlkCMMcbkiWsJREQiRGSeiGx1PsOzKfeiiKwXkY0i8rqIiLO/tYisE5FtmfcbY4wpHG72QJ4AFqhqPWCBs/0rInI10BFoBjQB2gBdnMP/Ae4F6jlLn0KI2RhjjMPNBNIPmOSsTwL6Z1FGgSAgAAgE/IHDIlIZKKuqP6vnQZYPs6lvjDHGS9xMINGqehDA+Yy6tICq/gx8Bxx0ljmquhGoCuzLVHSfs+83ROReEYkTkbgjR44UcBOMMab08upsvCIyH6iUxaEROaxfF2gIVHN2zRORzsC5LIpn+Ui9qr4LvOuc74iI7M7Jd2ehInA0j3WLmpLSlpLSDrC2FFUlpS35bUfNrHZ6NYGoao/sjonIYRGprKoHnSGphCyK3QQsUdUzTp1vgfbAR/ySVHDWD+QgnsjcxH9JvHFZPcpfHJWUtpSUdoC1pagqKW3xVjvcHMKaCQx21gcDM7IoswfoIiJ+IuKP5wL6RmfI67SItHfuvhqUTX1jjDFe4mYCGQ30FJGtQE9nGxGJFZEJTpmpwHZgHbAGWKOqXznHHgAmANucMt8WYuzGGFPqufZGQlVNBLpnsT8OuMdZTwPuy6Z+HJ5bewvLu4X4Xd5WUtpSUtoB1paiqqS0xSvtKFXTuRtjjCk4NpWJMcaYPLEEYowxJk8sgeSCiIwSkbUislpE5opIFbdjyisRGSsim5z2TBeR8m7HlBcicqszV1q6iBTL2y1FpI+IbHbmdfvNlD7FhYi8LyIJIhLvdiz5ISLVReQ7Z/699SLysNsx5ZWIBInIMhFZ47Tl6QI9v10DyTkRKauqp5z1h4BGqnq/y2HliYj0AhaqaqqIjAFQ1WEuh5VrItIQSAfeAf7u3FxRbIiIL7AFz52I+4DlwB2qusHVwPLAecj3DPChqhbmDS4FynkurbKqrhSRMGAF0L+Y/m8iQKiqnnEehfgf8LCqLimI81sPJBcykocjlGyefi8OVHWuqqY6m0v49YOZxYaqblTVzW7HkQ9tgW2qukNVLwBT8MwTV+yo6iLgmNtx5JeqHlTVlc76aSBj+qRiRz3OOJv+zlJgv1uWQHJJRJ4Tkb3AncBIt+MpIH/EnqNxS1Vgb6btbOd1M4VPRGKAlsBSdyPJOxHxFZHVeGb7mKeqBdYWSyCXEJH5IhKfxdIPQFVHqGp1YDLwoLvRXt6V2uKUGQGk4mlPkZSTdhRjWb3Hptj2bEsSESkDTAMeuWT0oVhR1TRVbYFnlKGtiBTY8KJrDxIWVZebv+sSHwPfAP/yYjj5cqW2iMhg4Hqguxbhi2G5+N+kONoHVM+0naN53Yx3OdcLpgGTVfULt+MpCKp6QkS+x/PupAK50cF6ILkgIvUybd4IbHIrlvwSkT7AMOBGVU1yO55SbDlQT0RqiUgAcDueeeKMS5wLz+/hmXfvFbfjyQ8Ricy4w1JEgoEeFODvlt2FlQsiMg24Cs9dP7uB+1V1v7tR5Y2IbMPzkq5EZ9eS4nhHmYjcBLwBRAIngNWq2tvdqHJHRPoCrwK+wPuq+pzLIeWJiHwC/A7P1OGHgX+p6nuuBpUHItIJ+BHPHHzpzu7hqjrLvajyRkSa4Xlhny+eDsNnqvpMgZ3fEogxxpi8sCEsY4wxeWIJxBhjTJ5YAjHGGJMnlkCMMcbkiSUQY4wxeWIJxJgCJCJnrlzqsvWnikhtZ72MiLwjItudmVQXiUg7EQlw1u1BYOMqSyDGFBEi0hjwVdUdzq4JeCYnrKeqjYE/ABWdSRcXALe5EqgxDksgxniBeIx15uxaJyK3Oft9RORtp0fxtYjMEpFbnGp3AjOccnWAdsA/VDUdwJmx9xun7JdOeWNcY11gY7zjZqAF0BzPk9nLRWQR0BGIAZoCUXimCn/fqdMR+MRZb4znqfq0bM4fD7TxSuTG5JD1QIzxjk7AJ85MqIeBH/D84HcCPlfVdFU9BHyXqU5l4EhOTu4klgvOC4+McYUlEGO8I6tp2i+3H+AcEOSsrweai8jl/o0GAsl5iM2YAmEJxBjvWATc5rzMJxLoDCzD80rRAc61kGg8kw9m2AjUBVDV7UAc8LQzOywiUi/jHSgiUgE4oqophdUgYy5lCcQY75gOrAXWAAuBx50hq2l43gESj+c97kuBk06db/h1QrkHqARsE5F1wHh+eVdIV6DYzQ5rShabjdeYQiYiZVT1jNOLWAZ0VNVDzvsavnO2s7t4nnGOL4Ani/n74E0xZ3dhGVP4vnZe8hMAjHJ6JqjqORH5F553ou/JrrLz4qkvLXkYt1kPxBhjTJ7YNRBjjDF5YgnEGGNMnlgCMcYYkyeWQIwxxuSJJRBjjDF58v+BAIHxSlzdpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot CV误差曲线\n",
    "test_means = grid2.cv_results_[ 'mean_test_score' ]\n",
    "test_stds = grid2.cv_results_[ 'std_test_score' ]\n",
    "train_means = grid2.cv_results_[ 'mean_train_score' ]\n",
    "train_stds = grid2.cv_results_[ 'std_train_score' ]\n",
    "\n",
    "# plot results\n",
    "n_Cs = len(Cs)\n",
    "number_penaltys = len(penaltys)\n",
    "test_scores = np.array(test_means).reshape(n_Cs,number_penaltys)\n",
    "train_scores = np.array(train_means).reshape(n_Cs,number_penaltys)\n",
    "test_stds = np.array(test_stds).reshape(n_Cs,number_penaltys)\n",
    "train_stds = np.array(train_stds).reshape(n_Cs,number_penaltys)\n",
    "\n",
    "x_axis = np.log10(Cs)\n",
    "for i, value in enumerate(penaltys):\n",
    "    #pyplot.plot(log(Cs), test_scores[i], label= 'penalty:'   + str(value))\n",
    "    plt.errorbar(x_axis, -test_scores[:,i], yerr=test_stds[:,i] ,label = penaltys[i] +' Test')\n",
    "    plt.errorbar(x_axis, -train_scores[:,i], yerr=train_stds[:,i] ,label = penaltys[i] +' Train')\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel( 'log(C)' )                                                                                                      \n",
    "plt.ylabel( 'logloss' )\n",
    "plt.savefig('Logisticgrid2SearchCV_C.png' )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3  使用LogisticRegressionCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegressionCV(Cs=[0.001, 0.01, 0.1, 1, 10, 100, 1000], class_weight=None,\n",
       "                     cv=5, dual=False, fit_intercept=True,\n",
       "                     intercept_scaling=1.0, l1_ratios=None, max_iter=100,\n",
       "                     multi_class='ovr', n_jobs=None, penalty='l1',\n",
       "                     random_state=None, refit=True, scoring='neg_log_loss',\n",
       "                     solver='liblinear', tol=0.0001, verbose=0)"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV\n",
    "lrcv_L1 = LogisticRegressionCV(Cs=Cs, cv = 5, scoring='neg_log_loss', penalty='l1', solver='liblinear', multi_class='ovr')\n",
    "lrcv_L1.fit(X_train, y_train) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: array([[-0.69314718, -0.6359066 , -0.48921992, -0.48815795, -0.48904032,\n",
       "         -0.48914202, -0.48915269],\n",
       "        [-0.69314718, -0.64460736, -0.52487947, -0.52917478, -0.53103422,\n",
       "         -0.53124003, -0.531256  ],\n",
       "        [-0.69314718, -0.63328109, -0.45825304, -0.45561705, -0.4558958 ,\n",
       "         -0.45592723, -0.45592845],\n",
       "        [-0.69314718, -0.63024185, -0.43335969, -0.42237952, -0.42195909,\n",
       "         -0.42192622, -0.42191763],\n",
       "        [-0.69314718, -0.63580105, -0.48687649, -0.48452611, -0.48410588,\n",
       "         -0.48411169, -0.48410863]])}"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lrcv_L1.scores_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
